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ABSTRACT 

Digital computer simulation has become an important 
technique for the evaluation of a wide range of systems. 

A general method, suitable for use in digital computer 
simulations, is often needed to measure the performance of 
those systems. A useful criterion for performance evalua- 
tion is the mean-square error between waveforms at two 
points in a system or between the waveform at one point and 
a theoretical reference. Signal-to-noise ratio, a useful 
parameter in many cases, can be defined for arbitrary wave- 
forms using measurements of power and mean-square error. 

Other system parM»ters, such as gain and delay, can be 
given general definitions through the use of a minimum 
mean-square error criterion. Techniques are developed to 
estimate power, gain, delay, slgnal-to-nolse ratio, and 
mean-square error in digital computer simulations of low- 
pass and bandpass systems. The techniques are applied to 
analog and digital communications. The signal-to-noise 
ratio estimates are shown to be maximum likelihood esti- 
mates in additive white Gaussian noise. The methods are 
seen to be especially useful for digital communication 
systems where the mapping from the signal-to-noise ratio 
to the error probability can be obtained. Simulation re- 
sults show the techniques developed to be accurate and quite 
versatile in evaluating the performance of many systems 
through digital computer simulation. 
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I . INTRODUCTION 

pigital computer simulation is becoming an increasing- 
ly important tool in the analysis of many systems. In 
order to evaluate systems through the use of computer simu- 
lation, a general quantitative measurement of performance 
that can be applied to a broad range of systems is needed. 

In many cases, an appropriate measurement of performance is 
the mean-square error between waveforms which appear at 
different points within a system, or between the waveform 
at one point and some theoretical reference waveform. In 
cuBBS irtiere the mean-square error itself is not an appro- 
priate figure of o»rit, other useful measures of performance 
may often be derived from the mean-square error. 

One qttantlty commonly desired in evaluating many sys- 
tems is the signal-to-noise ratio. By defining noise as the 
mean-square error, it is possible to use measurements of 
power and mean-square error to obtain signal-to-noise ratios 
for arbitrary waveforms. This procedure is complicated by 
the fact that many systems contain unknown gain and delays 
which need to be estimated in order to calculate the mean- 
square error. These parameters, also, can be defined through 
determination of estimates which minimize the mean-square 
error or maximize the signal-to-noise ratio. One area in 
which measurements of power, gain, delay, and the signal- 
to-noise ratio are particularly needed is in the evaluation 
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of analog and digital communicat ion systems through digital 
simulat it>n. 

The widespread use of eomplex digital eoinmim i ra 1 i on 
systems, and the advanrenu'nt I't eomimtt'r teehnology in re- 
eent years has resulted in t lu' inetu'asi'd use of digital eent- 
puter simulat ion t\>r prtnlieting t h«' per t ormaiiet' of digilal 
Oi'nimun i ca t ion systt'ms. Mi'iite t'urli' t t'l’hn i ipu's ha\'e bt'iMi 
applied to a wide* variety *'t syst»'ms, Obst'rvat icuis c>t t ht' 
simulated system under various operating conditions are used 
to predict the performance of the actual system. 

The most common method for evaluating the performance 
of a digital conoaunicatlon system through computer simula- 
tion is to determine the symbol error probability by direct 
error counting. However, this technique has the disadvantage 
of requiring prohibitively large amounts of computing time 
to measure typically encountered error probabilities. For 
example, on the average, to produce ten errors at an error 
probability of requires that lo' information symbols 

be simi ated. Extremal statistics and rt'lafed tcchnique.s 
can generally extend the practical range of the error 
counting method by apprviximat el y an ordt'r of magnitude. 

However, estimation of error probabilities below about 
4 — S 

10~ or 10 ' still require too much computer time to be 
economically feasible in most cases. 

When the mapping from the signal-to-noi se ratio to 
the error probability is either analytically known t.r can 
be experimentally determined, an alternative approach is 


to estimate the signal-to~noise ratio and map this to the 
symbol error probability. In the work which follows, a 
method is developed to accurately measure the signal-to- 
noise ratio in a digital computer simulation of a communi- 
cation system. The method is quite general and may be 
applied to a wide variety of systems. Both digital and 
analog systems can be evaluated. The use of direct simula- 
tion is not restricted to systems with lowpass signals. 
Simulation of a bandpass system may utilize either conver- 
sion to an equivalent baseband system, in which lowpass 
signals are processed, or a direct simulation of the system, 
which retains signals having bandpass spectra. The latter 
approach is coninonly referred to as direct RF simulation, 
even when the signals are not radio frequency. 

In the work which follows, computer routines are de- 
veloped for implementing signal-to— noise ratio estimators 
for the lowpass case and for the bandpass case where direct 
RF simulation of a bandpass system has been utilized. The 
lowpass estimator forms a minimum mean-square error estimate 
of the signal, allowing for any unknown amplitude scaling 
and time delay within the system. In the bandpass case, 
both group delay and phase delay of the signal may be 
arbitrary. The FFT provides a practical method for esti- 
mating delay in the lowpass case. For direct simulation of 
bandpass systems the FFT is essential for independently 
determining group and phase delay to form the minimum mean- 
square signal estimate. An option is provided for applying 
a data window to smooth the estimates if desired. 
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To evaluate the accuracy, reliability, and versatility 
of the signal-to-noise ratio estimation routines, a variety 
of tests are performed. First simple filtered iu'>ise tt'sts 
are used to verify the accuracy and deterniiiu' t lu' useful 
range of the est iniati>r.s. Then digital ci'niputt'r .si mu 1 at i<.>n 
and performance eva luat i ».'ns are ci'ndui'tt'd for a numlH'r ot 
systems including FM systt'ms with pliast'- 1 oi'keii loi^p d^'iniHlu 1 a- 
tion. coherent ASK and FSK. noncoherent FSK, and RayUMgh 
fading channels, including diversity transmission. The 
simulations incorporate additive Gaussian noise, and signal- 
to-noise ratio estimates are obtained to beyond 50 dB in some 

3 

cases using simulations of approximately 10 information 
symbols. The simulations were obtained using the IBM 370 
computers at the University of Missouri and a Univac 1110 
conqputer at NASA-Johnson Space Center in Houston. Texas. 

Results obtained from the estimators are excellent. 

In all the various systems which were simulated, accurate 
results were obtained over the ranges of the s i gnal - 1 o-tu' i se 
ratio that were considered. Of course, the necessity I’f 
knowing the mapping from the signal-to~noise rat it' to the 
symbol error probability is sometimes a handicap in dealing 
with digital systems. However, in many cases this mapping 
can be determined using analysis or approximately deter- 
mined using bounding techniques. Sometimes, however, the 
signal-to-noise ratio is a satisfactory figure of merit, 
and determination of the symbol error probability is un- 


necessary . 
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The methods developed appear to be quite versatile 
and promise to provide a very useful tool for evaluating 
the performance of a wide variety of practical systems 
through the use of computer simulation. 
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II. REVIEW OF SIMULATION TECHNIQUES 

Computer simulation of communication systems is bt»- 
coming a popular technique in the study of systems which 
are too complicated for an exact mathematical analysis. 
Simulation offers an economical alternative to th(' fabri- 
cation and testing of all the variations of the system 
being studied. In order to take advantage of these com- 
puter methods, several schemes have been developed for 
evaluating the performance of systems through the use of 
eoogmter simaiation. 

Host c oCT Bi u nicat ion systems require some modulation 
technique which generates a bandpass signaling scheme. The 
most direct method of accurately simulating such systems 
is to generate simulations based upon these bandpass signals. 
This method has been used successfully by several investiga- 
tors [1] - [6]. An alternative approach is to mathe- 
matically convert the bandpass signals to equivalent low- 
pass (baseband) signals. The baseband technique avoids 
some problems encountered in directly simulating the band- 
pass systems and often lends itself to simpler analysis. 

For these reasons, the baseband methods have found somewhat 
wider acceptance than direct RF simulation of bandpass sys- 
tems and appear more often in the literature [7] - [12]. 


All numbers shown in brackets refer to corresponding num- 
bers in the list of references. 
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One advantage of baseband simulation Is that often 
less samples need to be used lor a given number of informa- 
tion symbols. This results from the fart that the highest 
frequencies being sampled are those In lowpass signals. 

For direct RF simulation of the bandpass signals, the 
sampling frequency must be sufficiently high to includi 
carrier frequency and usually an upper sideband of modula- 
tion. The number of samples necessary for a given Simula 
tlon time Interval is generally at least live to ten times 
less with the baseband method. Another advantage la that 
the Btmuletloa models are sometimes easier to program. For 
esample. a bank of filters written for use In an exclusively 
baseband simulation may often Include only lowpass filters, 
in the direct RF simulations, the usually more complex 
bandpass filters are aecessary . Problems such as adjaceat 
channel Interference can be eliminated in baseband methods. 


where no adjacent frequency channels exist. This can be an 
advantage if such an Idealized system is to be simulated. 

More often, however, this inability to directly account 
for interference is a disadvantage. Many times in an FSK 
system or in a frequency multiplexed system, the problem of 
Interference between adjacent frequency channels is an impor- 
tant effect which must be included in the simulation. The 
direct RF simulation can simply model the system as it 
actually exists, simulating all channels simultaneously at 
the appropriate carrier frequencies. The baseband simuli- 
tion is not capable of this approach, so more elaborate tech 

niques are required. 
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One simulation method used by Hedderly and Lundquist 
[13] to account for the interference effects of frequency 
diversity is to compute the spectrum of power from adjacent 
channels which would fall into the channel of interest and 
then to Increase the additive channel noise to include the 
additional "noise" power created by Interference. Another 
technique used by Bello and Ehrman [14] involves the use of 
short tone bursts of the desired frequency channel adjacent 
in time to tone bursts of two Interfering frequencies. The 
channel is used in rapid time diversity to account for the 
late»yi«»l interference effects In an PSK system. These 
methods have been applied In BF simulations with some 


success, but the direct modeling of all the frequency 
channels is certainly more intuitively appealing. 

Another disadvantage of the baseband technique is that 
it is usually necessary to model inherently bandpass systems 
using a complex envelope representation. Hence, either two 
real number simulation paths or a complex number simulation 
must be used to represent the general system. 

In the work which follows, the direct PF simulation 
approach for modeling bandpass systems was most often used. 
It was felt that the baseband simulations often Idealized 
and oversimplified some of the various effects which could 
not be directly modeled. Also, intuitive insight into a 
system seems to be more easily applied with direct RF model 
ing. As explained later, the accuracy of the method chosen 
to evaluate the performance of the systems being simulated 

total number of samples used and not 


is a function of the 
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of the number of modulation symbols involved. Hence tht» 
increased number t>t sample.s requirt'd by tht' direct RF simu- 
lation techniqut' ti> simuiatt' a ^^ivt'ii number of mtHiula(ii>n 

symbi>ls compart'd 1 1 » bast'baiui mt'tlu>ds is lU't a serious 

prt'b 1 t*m . 

Oiu'e a siimilat it>n. bast-baiui or KK. has bi'eii lomplt'ted 
and data ctillt'ctcd. s»)nu> nu'thod must bt> ust'il to I'valuate 
the system performance. In certain limited applications, 
only qualitative information is required. Simply plotting 
output waveforms and inspecting them visually for distortion 
is sofltetines acceptable in these cases [153. However, la 
most cases a quantitative measurement of system performance 
is required. For digital communication systems, the basic 
measure of performance is usually taken to be the probability 
that a given information symbol is in error. Several ap- 
proaches have bt'en applied to directly or indirectly «'stiniate 
this quantity. The most straightforward method i.s simply to 
periorm the simulation, count the errors, and utilize the law 
of large numbers to directly estimate the probability of 
error. Other methods use extremal statistics and related 
techniques to estimate the probability of error. Still 
another possibility is to measure some other quantity such 
as the signa 1 -to-noise ratio and estimate the probability 
of error indirectly. 

By far the most often used technique in the estimation 
(the probability of symbol error) for digital communi- 
cation system.s has been direct error counting. This 
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approach has been used by many researchers and is reported 
in the literature [161 - llH). This approach works very 
well where it is practical. However, actual systems typi- 
cally display error probal. i I 1 1 i es ranginp from to lo"^ 

or evt'n It'ss. To obtain a rt'liable t'stimatc ot t hi- ,'rror 
probability. it is desirable to utili.a> at least eiiou^ib 
data to includt' ten I'rrors. For a siimilation ot a system 

e.xhibiting an error probability of 10“‘’. this would repuir«> 

7 

10 symbols to be processed. This is completely impractical 
from an economic point of view. For example, Leon and 

Kitaharft [19} report that for a particular ^K ( quadrature 
friiaao shift keyed ) simulation performed on a OW: 6500 com- 
puter, the expected cost of generating a single error with 
an error probability of 10 ® would b<- approximately $4000. 
They obtained estimates down to about 10~^ with direct 
counting methods. Similar results were reported hy I.eim 
and others [20] for PSK (phase shift keyt'd) data links. 

Bello and Crystal [21 1 were able to use dirt'ct erretr 
counts for P^, to about 10~‘* . Their system was one which 
was difficult to analyze for low signa 1 -to-noise ratios. 

For high signal-to-noiso ratios, analysis was mathematically 

tractable and bounds were developed for P,, below 10~"^ 

F, 

Joining the bounds with the computer simulation results 
yielded curves over a wide range of Pj, . Because of the 
large amounts of computer time required, the direct error 
count method seems to be impractical below error probabil- 
ities of .about 10 in most cases, or lO”'' in the most 
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efficient siniulatlons. This Is not adequate for most appli- 
cat ions . 

in an effort to estimate low error probabilities using 
less data, several researchers have applied estremal statis- 
tics to the problem [221 - | 26 1 . m this technique an 
exponential function is used to approximate the tail of the 
density function for the errors by adjusting parameters in 
a general exponential form. The parameters are computed 
using the statistics of the maxima (extremes) of groups of 
independent samples collected from the simulation [ 27 ], 
Smaller values of Pj are then computed from this easily 

e.3K»«tl.l awro.lm.tloo to the density func- 
tion. This general technique may be applied to a wide 
variety of density functions which behave as a decaying 
exponential on the tails. Specific tests are available to 
determine whether a set of data points do in tact satisfy 
this type of behavior [28], One disadvantage of the ex- 
tremal statistics method is that it requires a learning 
period when the signal is known so that the parameters may 
be computed as a function of only the unknown noise dis- 
tribution. A method to eliminate the learning period and 
estimate parameters based upon noisy samples has been 
proposed by Milsteln [29]. Schwartz [30] and other in- 
vestigators have concluded that extremal statistics may be 
used to reduce by an order of magnitude in many cases the 
data that are required tor satisfactory Pj. estimates. 

Attempts to extend these methods beyond about two orders of 
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inugnitudo oft«?n rosult in Pj,, ostimatt^s which urc in error 
by almost an order of magnitude, howt'ver . It app**ars 
that extremal statistics may t'xtend t lu‘ capabilities of 
t ht' t'ri'i'r i'i>unt i ng t t'l’hn i ipu's I rt>m a typii’al limit 

to about 10 . or at best 10 * . This may be sutticient 
in si'iiu' a pp I i »-■ a t i ons . but it is still tar short ot the 
or 10~^ range which is oftt'ii dt'siriul. 

Another variation of error counting has been tried 
with some success by a number of researchers. This method, 
described by Gooding [31], is based upon measuring what is 
temed a ’’pseudoerror'* probability and mapping this quantity 
to Pp. To implement this method, thresholds are established 
for the demodulated signals other than the usual decision 
thresholds used to decode the received signal. These new 
thresholds define a region similar to that which corresponds 
to an erasure in the classical binary erasure channel. All 
demodulated symbols whic*h full inti> this highly uncertain 
region within the now thresholds are termed pseudoerrors . 

The symbols may still be decoded as usual; thert> is no need 
to actually erase them. By monitoring tht'st* pseudiu>rrors 
and applying knowledge i>f the distribut ii>n luui’t ii>n whic’h 
generates them, Pj. may bt' estimated. A signiticuiit advan- 
tage of this technique is that it requirt's no knowledge of 
the input data stream. The thresholds are gonorally estab- 
lished to be symmetric with respect to the ideal signal 
points in decision space, so that the pseudoerror probabil- 
ity is independent of the stream of input data and even of 
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its statistics. This method does, however, require that the 
distribution function for the noise be known in order that 
proper thresholds may be established and the mapping to the 
error probability computed. The method has been applied tt> 
PSK and QPSK systems [ 32 I - [3^1. It appears t«> yield re- 
sults slightly more reliable than those obtained by extrc'mal 
statistics methods. However, it should be remembered that 
the extremal statistics approach does not require that the 
distribution function be precisely known. The error prob- 
ability estimates obtained using the pseudoerror technique 
are usually accurate to within a factor of three for values 
of P-c as small as two orders of magnitude less than those 

which are practical with direct error counting. This method 

-5 

is not usually capable of measuring Pg much below 10 in 
computer simulation work, and Leon and Kitahara [35] note 
that this technique is not economically feasible for what 

_7 

they term practical 10 error probabilities. 

Another similar technique has been suggested by 
Weinstein [36], which involves counting the decoded symbols 
which fall above some threshold. This method also is 
reported to be capable of reliably estimating Pg down to 
values between one and two orders of magnitude smaller 
than is possible with a direct error count. The technique 
is based on linearizing the tail of the distribution func- 
tion of the noise with a log-log transformation. 

An indirect method for estimating the error probability 
is to first estimate the signal-to-noise ratio and map this 
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into Pp. This method has the disadvantage that the dis- 

£1 

tribution function for the noise must be known, or somehow 
the mapping from the signal-to-noise ratio to the prob- 
ability of error must be obtained. The estimation of Pj,, 
through estimation of the signal-to-noiso ratio has ont' 
striking advantage over all the other techniciues mentiont'd. 
For any signal-to-noise ratio that can bt' reliably estimated 
where the mapping to Pg is known, reliable estimates lor the 
error probability can be obtained. That is, there is no 
reason why P*. cannot be measured for any range of interest. 
For most digital systems, if signal-to-noise ratios as high 
ns 30 or 40 dB can be measured, probability of error estim- 

-7 

ates are possible far below the typical 10 range. 

A B^thod used for estimating the signal-to-noise ratio 
fB order to monitor the error rate of a PCI! syst^ is given 
by Gagliardl and Thomas [371. The method presented is only 
applicable to a system or a simulation with no unknown delay, 
however. No implementation is suggested for use with a 
computer simulation. The decoded data stream at the output 
is assumed to be correct in order to serve as a reference. 

If the method is to be applied exclusively to systems with 
very low error probabilities, this might be acceptable. It 
^ifould seem more general and more accurate to use a known 
input data stream if possible. This does, however, intro- 
duce the problem of finding an unknown delay between input 
and output, a problem not treated in [38]. 
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Another method for signal-to-noise ratio estimation is 
suggested by Nahi and Gagliardi [39]. For the case where 
both signal and noise are Gaussian processes with non- 
identical correlation functions known to within a constant, 
the signal-to-noise ratio may bo estimated using a hard- 
limiter. The method is based on known properties of the 
correlation function of the output of a hard-limiter. This 
estimator is clearly very specialized, and its reliability 
has not been analyzed. It is not a suitable estimator in 
general and is mentioned only for cooqpleteness. 

the signal>to>noise ratio estimation technique for find- 
ing Pg would se^ to be the most promising method for small 
error probabilities, if some method is available for ac- 
curately measuring the signal-to-noise ratio within a digital 

-5 

simulation. Indeed for Pg below 10 , it seems to be the 

only practical method, since all other techniques are not 
economically feasible because of the large amounts of com- 
puter time they require. It has been suggested [40] that 
no single technique for Pg estimation should be intended 
for a very wide dynamic range of error probabilities and 
that different methods should be applied on different ranges. 
Whether this is true or not, the missing range for Pg estim- 
ation through simulation seems to be exactly the range most 
needed for computer simulation of practical systems, the 

_5 

range where the error probability is below 10 

The method used in this work for predicting the per- 
formance of digital communication systems through computer 


16 


simulation is estimation of the signal-to-noise ratio, 
which is mapped to the error probability. Thus the map- 
ping from the signal-to-noise ratio to must be known, 
or at least approximated. Also, a representation of the 
signal must be known to serve as a reference. However, 
it is assumed that the system may Incorporate any unknown 
amplitude scaling and any unknown linear phase characteris- 
tic in processing the signal, so that the exact form of the 
signal at the point of measurement may be quite different 
from the form of the reference. None of the other techni- 
ques descrtbed shove require knowledge of both the slgasl 
form and the mapping from the slgnal-to«nolse ratio to 
^wever, each method requires that either the signal or the 
noise distribution function be known. The big advantage of 
the signal-to-noise ratio method Is that: practical systems 
operating with practical error probabilities can be analyzed 
through computer simulations. No other technique is capable 
of this. And it is worthy of note that in this process 
for estimating Pg, an estimate for the signal-to-noise 
ratio is always obtained. The signal-to-noise ratio 
itself is often a useful parameter in evaluating communi- 
cation systems without reference to an error probability. 
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III. DESCRIPTION OF THE PROBLEM 

If the SNR (signal-to-noise ratio) is to be adopted 
as a parameter for evaluating the performance of all sys- 
tems to be analyzed, a gt'm'ral definition for this quant ily 
must be developed. For most practical (.‘(.)mmun i I'at ion sys- 
tems, the amplitude response and phase <'haracter i st . c ot tht' 
system are not precisely known. Therefore, at any arbitrary 
point in a system where an estimate for the SNR is desired, 
it is often difficult to completely specify the signal. 

Tbv^r Itefore tbe SHE can be estlBated, a rigorous definitiem 
for tbe signal must be devised. 

It is usually possible to provide an input to a system 
so that the signal coB^)onent of the waveform at the test 
point is known except for changes imposed within the sys- 
tem by an unknown scale factor and linear phase characteris- 
tic. A linear system is usually considered to be distortion- 
less if it imposes only amplitude scaling and a time delay 
upon an input signal [41]. It seems appropriate then to 
allow for any arbitrary scale factor and delay to be present 
in the signal component at the point where the SNR is to be 
estimated. For lowpass systems, this signal model is 
adequate . 

In such a case, Shepertycki [42] has suggested an 
error measurement scheme for telemetry systems. In his 
method of error measurement, the mean-square error is com- 
puted between the waveform to be analyzed and a reference 
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signal. The amplitude scaling and time delay of the ref- 
erence are adjusted to minimize the mean-sQuare error, and 
this error is used as a figure of merit for the system. 
Simpson and Houts [43] have proposed that the delay value 
which minimizes the mean-square error be defined as the 
average time delay for a linear sy.stem. Later, the .same 
investigators [44] suggested that the minimum mean-.square 
error technique can be used in the analysis of waveform 
distortion in linear systems. It is also suggested that 
random noise could be included in this analysis. 

When a bandpass syst^ is being analyzed using direct 
RF siOfUlatiOB, both group delay and phase delay are in- 
volved. A more general definition for the signal coo^onent 
is necessary for this more complicated case. In bandpass 
systems involving modulation, the Information is carried in 
the complex envelope of the signal. To avoid corruption 
of the information, it is sufficient that this envelope 
remain undistorted^ delay of the carrier itself is not 
important in preserving the information content of the sig- 
nal. It is possible to model phase delay as a time delay of 
the carrier, and group delay as a time delay of the complex 
envelope. Therefore, a definition of the signal in this 
case should allow for any arbitrary scale factor, delay of 
the envelope (group delay), and delay of the carrier (phase 
delay) to be present at the point of measurement. 

In keeping with this reasoning, a definition was de- 
veloped for the signal component of any arbitrary bandpass 
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waveform to be analyzed. It is assumed that a reference 
signal is available which differs from the signal com- 
ponent only by an unknown amplitude scale lai’tor and un- 
known group and phase delays. P'or lowpass systi^ms th(' 
group and phase delays are equal, and only the scal(' fac- 
tor and a simple time shift are unknown. To del'int' the 
signal at the point of measurement, the mean-square error 
between the reference signal and the waveform to be analyzed 
is observed. The form of the reference after amplitude 
scaling and appropriate time delays have been imposed upon 
it to miniflilae the mean-square error is defined to be the 
signal conq>onent of the waveform under test . All other 
components of the waveform are grouped together and are 
termed noise. Thus errors resulting from both random noise 
and system distortion terms are included. With this defini- 
tion, it is possible to define the SNR for any system where 
an appropriate reference signal can be obtained. In digital 
computer simulation, such a reference waveform is almost 
always available. 

Once the reference and noisy waveforms have been ob- 
tained, the problem of estimating the SNR involves finding 
a method to determine the gain and delay values which yield 
the minimum mean-square error. It will be shown in later 
sections that the appropriate gain and delay values can be 
obtained using the cross-correlation function between the 
noisy waveform and the reference. For bandpass systems, the 
appropriate value of group delay will be seen to be the value 


f vrifil .. .. ^T r*. i ; ■::■ u.d. 
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for which the envelope of the oross-eorrelat ii>n tunetion 
is maximized. For a lowpass system tite peak of the itoss- 
oorrelation fmut it>n itstMf iiuiieates t ht« averag,' system 
delay. In either ease t ht' eross-ei>rre I a t ion liuutii>n must 
he found. Fi>r lurpe sets of data, the eomputt-r time re- 
quired to direetly K«*tO‘rate this funetion is prohibitive. 

By utilizing the iFT, the eross-eorrt' I a t ion luru-tion «’an 
be computed quickly. Hilbert transform methc^ds provide a 
way to obtain the envelope of thi.s function when it is 
required. Scanning these functions for peak values provides 
a delay estimate. However, it was experimentally determined 
that the errors introduced in the successive transforma- 
tions BMke the values obtained for the cross-correlation 
function unacceptable for estimating the SNR. except in the 
very noisy cases where the SNR is quite low. To obtain a 
good SNR estimate, a more accurate method is required for 
the final computations. 

In the lowpass case, the signal component of the wuvt'- 
form may be obtained by merely time shifting and amplitudt* 
scaling the reference signal. A relatively simpli> method 
was found for accurately estimating the SNK in this case. 
When the system is bandpass, arbitrary group and phase 
delays are both involved; simple time shifting will not 
account for these effects. No simple scheme was found for 
deriving the signal component by performing a trivial opera- 
tion on the reference for this case. The most direct 
methods seem to involve either a ct>mpl i la t,'d filtering 
operatii>n or .suci'ess i vt' Fourit*r t runs format i i>ns . Tlu' 
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errors introduced in these methods defeat the purpose, which 
is to obtain a highly accurate representation of the signal. 

A method was found to measure the SNR in direct RF simulations 
of bandpass systems, but it is more complex than the lowpass 
case. Therefore, lowpass and bandpass systems are treated 
separately in the following development of the estimators. 

A question that is considered is that of how accurate 
and reliable the estimates obtained for the SNR are. A 
partial answer is obtained by simply using the estimators on 
a variety of systems where the SNR can be estimated by 
another means for comparison. Consideration is also given 
to the statistical properties of the estimators. It is help- 
ful to compare these estimators to maximum likelihood estim- 
ators. Confidence intervals to indicate the reliability of 
the estimates are developed where practical, also. 

After the SNR estimators have been developed and imple- 
mented, they are tested on a wide variety of systems. Both 
analog and digital systems are simulated and evaluated using 
the SNR estimators. The systems are necessarily chosen to be 
suitable for theoretical analysis, since the estimators 
themselves are being evaluated. For digital systems the 
desired parameter for evaluating performance is usually the 
symbol error probability. Therefore, in the testing pro- 
grams, digital systems are simulated where the mapping from 
the SNR to Pg is known. For Pg ranges where it is prac- 
tical , a logical test for digital systems is to compare the 
error probability estimate computed from the measured SNR 
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to the error probability obtained by direct error count. 
Obviously this testing is restricted to large error prob- 
abilities, since the direct counting technique' imprac- 
tical elsewhere. However, it provides an e.xcellent test 
for the SNR estimators where it can be applied. Where the 
SNR is large, it is only possible to directly test the 
estimation of the SNR itself, using simple systems where 
gain and phase information are known, and the true SNR can 
be computed. For this range it is not practical to make 
direct c<»iparisons with the error probability, since there 
exists BO feasible alternative to the SNR estimators for 
determining Pg at low error probabilities — the very reason 
for developing these estimators in the first place. 
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IV. AN SNR ESTIMATOR FOR 
COMPUTER SIMULATIONS OF I,0\VPASS SYSTEMS 

^ • THEORETICAL DEVELOPMENT 

In nnuly/in^ a o«.)n*nun i cat ion systi'm by computt'r Simula 
tion. or in ova lua ting t lu' simulation it sol I, quant itios of 
intorost aro t lu' moan-scjuaro terror and t ho s i una 1 - 1 o-nt' i so 
ratio at various points in tho simulation. Measur«>monts of 
system gain, system delay, and power are also important 
in the analysis of many systems. Subroutine SNRMSE is a 
computer routine which was developed to estimate these quan 
titles In digital computer simulations of wide-sense sta- 
tionary lowpass systems. 

In computing all these quantities except power, two 
sets of data are required: a reference array and an array 

of measurement data to be analyzed. For the computation 
of the error, the measurement data is compared to an 
amplitude scaled and time-shifted representation of the 
reference data. A lowpass system which exhibits only- 
amplitude scaling and a time delay is considered to be 
distortionless or ideal. Optimal estimatt's for system 
gain and system delay are defined to be t ht>se values of 
the gain and delay which minimize the mean-square error 
between the measurement data and the reference data. 

The reference signal will be denoted by x(t). and the 
measurement data by y(t). Since the ideal system is 
assumed to impose only amplitude scaling and a time shift 
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on the reference, the ideal signal at any point may be 
represented by Ax(t-T), where A is the assumed gain of the 
system and t is the assumed delay. Then, letting an over- 

bar denote a time average, the mean-square error is givi>n 
by 

E(A.t) = ly(t) - Ax(t-i ) (1) 

Figure 1 shows the test configuration which is assumed for 
computing the mean-square error. Expanding (1), 

“ y^(t) - 2Ax(t>T) y(t) + A*K*<t-T> 

E(A,t) - y^(t) - 2Ax(t-x) y(t) + A^x^(t-T). 

And, assuming a stationary system, 

E(A,t) . Py - 2AR^y(x) t (2) 

where 


Py = power in the measurement data y(t) 

= power in the reference data x(t) 

~ cross-correlation function between the ref- 
erence data and the measurement data. 

It is clear that and A^P^ are positive quantities. 
Hence, for any value of A, E(A,t) is minimized when 

is maximized. Assuming system gain A is posi- 
tive, the optimal estimate for system delay must therefore 
be the value of delay for which maximized. 





mK, 
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Denoting this value of delay by the mean-square error 
becomes 


E(A.t,„) = P, - 2AR (f ) + A“P . 
m y xy' m' x 

fixed, the expression for E(A,t^^) is diflt'ren- 
tiated with respect to A and set equal to zero to yield 
the optimal estimate for system gain. This estimate is 
denoted by A so that 


<1E(A,t^) 

dA 


A-A. 




m 


+ 2A P » 0. 

ID X 


Thus 


A = 

m P 


(3) 


From (2) and (3), the resulting mean-square error is given 
by 


m m 


= P - 


(t ) 
xy m^ 


(4) 


From the definitions of A^ and the estimate of the 
signal at the point where the measurement data is taken 
must be Aj^x(t-r^). Thus the signal power at that point 
in the system is estimated to be 




A^P 
m X 


( 5 ) 
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or 


S = 





( 6 ) 


In obtaining the SNR estimate, all components ot y(t) 
which are not included in the signal estimate are termed 
noise. That is, the SNR is defined as the ratio of th<? 
signal power to the noise power at the point where the 
measurement data is taken; and noise will include any dis- 
tortion or roundoff errors present in addition to the 
additive random noise. Then the SNR may be expressed by 



or 


SNR = 


2 

R (t ) 
xy ' m'^ 

2 


P^P - 
X y 




< 7 ) 


If the normalized cross-correlation is defined by 


P 


x(t-T)y(tT 
A^(t) • /y2(t) 


R (t ) 

/P P 
X y 


then -15 p 51. It is possible to express the SNR in the form 


SNR 



1 


( 8 ) 
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(9) 


( 10 ) 


Estimates of the various system parameters in terms of 

Pjj» Py. given in Table I. 

Equation (8) illustrates that the SNR estimate is 

accurate if the aonnalized cross-correlation p is accurately 

determined. It is is^>ortaat to note that p is independent 

of the power contained in the reference and measurement 

data. Therefore, very accurate estimates can be obtained 

for the SNR, even in the presence of large errors in the 

individual measurements for P , P , and R (t). As long 

A y Ay 

as the same data are used to compute P and P as were 

A y 

used in determining R (x), errors in the measurements 

xy 

often tend to be cancelled when R (r) is divided by 

xy 

✓P P . This effect is used to advantage in the SNR 
X y 

estimation routine. 

When the SNR is to be estimated in a digital computer 
simulation, the values of Rjjy^^m^’ ^y 

estimated from a finite set of sample values. If a set 
of N samples are to be used, equation (7) becomes 
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TABLE I. 


PARAMETER ESTIMATES FOR LOMPASS SYSTEMS 
IN TERMS OF P^, Py . 
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(SNR) = 


H Ib <h> - M ■ 

^1 — X J^X”! 


( 11 ) 

where the subscript denotes a time sample. 

It is possible to determine whether this SNR 
estimate is a maximum likelihood estimate for some addi- 
tive noise channels. Consider the case where the mea- 
surement data are given by 


y(t) - A^x(t-Tj + n (t) 


( 12 ) 


The estimates described earlier for gain and delay will 
alwa 3 ^ give a representation of y(t) in this form where 
n (t) represents all components of y(t) which are ortho- 
gonal to x(t-Tjjj). This may be seen by noting that from 
equations (1) and (12) 


2 ^ 

n (t) = [y(t) - A^x(t-T)] = E(A^, t^) 


(13) 


and also from (12), 

- 2 j 

y^(t) = tAjjjX(t-Xjjj)]^+n' (t) + 2Ajj^x( t-tjjj)n (t). 
Substituting (5) and (13) into the above expression yields 


P * S + E(A_,t ) +• 2A x(t-T ) n (t). 
y m m m m 


( 14 ) 
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However, combining (4) and (6) shows that 

Py = S 4 E(A„.T^). 

Thus, comparing (1-i) and (15). 

In fact, an alternative derivation of the SNR estimator 
can be easily accomplished by using the orthogonality 
principle applied to linear minimum mean-square estimation 
Thus, the representation of y(t) in (12) is valid where 
n’(t) represents all components of y(t) orthogonal to 
x(t-T ). This representation is useful in determining 
whether the estimate given by (11). based on minimizing 
the mean-square error, is also the maximum likelihood 

(ML) estimate for the SNR. 

To find the ML estimate for the SNR it is first 

necessary to obtain the joint probability density func- 
tion for samples from y(t) as a function of gain and 
delay. For the additive white Gaussian noise channel 


y(t.) = Ax(t.-x) + n(tp < ^ 

where the subscript i denotes a time sample. The noise 
samples are independent and zero mean. Thus the samplts 
of y(t), denoted y(t^. are Gaussian with means given by 
Ax(tj-T). These samples are also independent. Denoting 
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the density function for y(t^) by f^Cy^^jA.i), the joint 
probability density function for N samples from y(t) is 
given by 


N 

f_(y |A,i) = n f.(y^|A,i) 

y 1 i=l I 


(17) 


The values of A and t for which this likelihood function 

is maximized are the ML estimates for gain and delay. If 

2 

the noise variance is denoted by o^, the likelihood func- 
tion becomes 


N 1 r / 

f (yU.T) - n exp - 

n 


y(t^Ax(tj-T) 

o 


or 


2 2 

f_ (y|A,T)»(2TT0^) exp^- 


1 N 

I 

2al i=l 


y^(t^) + A2x2(t.-x) 


- 2Ax(t^-r)y(tj) 


The log-likelihood function is 


N 


- — ^ [ [y^(t.) + A^x^Ct.-T) - 2Ax(t ,-T )y(t . )1 


20“ i=l 
n 


- I ln(2TTo^) 


or 






f 
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1 



N 2 

2A ^ x(t .-T )y(t - )-A 
i=l 


N , 

I X- 
i=l 


(t^-T )- 


N 2 

I y^tj) 

i = l ' 


- I ln(2.o2). 


Since x(t) is assumed to be wide-sense stationary, the 

second summation term is not a function of x. Thus it 

is clear that, assuming positive gain, the log-likelihood 

function is maximized ar a function of x when the term 
N 

J x(t.— T)y(t.) is maximized. From (7) and (11), the 
1«1 ^ ^ 1 ** 

discrete form for R^„(t) is ^ • X x(t .-x)y(t . ) . Therefore, 
the delay value required to maximize the log- likelihood 
function is Hence the delay estiMite given by the 

minimum mean-square error technique is the ML delay esti- 
mate for the additive white Gaussian noise channel. 

It can be shown that the ML estimate of the ratio 
of two parameters is simply the ratio of the ML estimates 
[45]. Therefore, in order to estimate the SNR it is only 
necessary to form ML estimates for the signal power and 
the noise power and confute the ratio. Denoting the 
power in the reference by P and noting that the discrete 
form for in (11) is ^ log- likelihood 

function may be expressed as 


I 

f 
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1 



N 


2A/F' 1 

^ i=l 

1 / 


y(t.) - NA^P^ 


N 

I 


i=l 


y^(t.) 


N 9 

- I ln(27iap. 

The power contained in the signal component of y(t) is 

given by A P^. Defining = aS^ and rewriting the log- 
likelihood function yields 

1 


- I 



^IVV' 


y(tj) - HS, - 


K 

I 

i-1 


y^(t. ) 


again assuming that A is positive. 

The ML estimates for and (o^) can be found by 

setting the partial derivatives with respect to S and 
2 ^ 

(a^j) equal to zero and solving. Then the ML estimate for 

the SNR is the ratio of these estimates. 

For the signal power, the partial derivative of the 
log-likelihood function yields 


1 N 
— I 




y(tj) 


- N = 0 
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Sx'?- 

X 


N Jj 


(18) 


For the noise power, (o^). the partial derivative of 


the log-likelihood function yields 


1 


2(0^)^ 


pr N 

;/s y 

X L- 


'X(t.-T ) 

1 tn'^ 


i=l /F 


y(t. ) 


N 


- NS^ - I 


i=l 




2(a2) 


= 0 


and ueiitg (18) 


[2NS 

(0=) * 


- NS - 


I y^t )] 

i»l ^ 


-N 


or 


1 ^ o 

% ’ I y <‘i) - 


(19) 


Then, expressing (SNR) as the ratio S /a^ yields 

X tl 


(SNR) = 


P 


fl ? 1' 

Kill ' 

5 ! y\t,) . S 


Using (18) to replace and using the sunraation which 
represents the discrete formulation for P results in 

X 


(SNR) = 




( 20 ) 


1?2 
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which is identical to (11). Thus the SNR estimator de- 
veloped using the minimum mean-square error criterion 

yields the ML estimate of the SNR for the white additive 
Gaussian noiso channel. 

If confidence intervals are to be found, it is n<‘c- 
essary to determine the probability density function of 
the SNR estimate. In the additive white Gaussian noise 
case, the required density function can be found unless 
the SNR is low. Except for low values of the SNR, it has 
been experimentally determined that delay is estimated in 
almost all cases with sufficient accuracy that the error in 
the SNR estimate generated by inaccurately determining 
delay is small compared to other error sources. When delay 
errors of significant size do occur, they are usually of 
sufficient magnitude to be quite obvious. Therefore, it 
IS assumed that x, the system delay, has been correctly 
determined. In order to determine the density function 

for the SNR estimate, first denote the term in brackets 
in (18) as 


1 N 

G = N I x(t -T)y(t. ). 

1=1 ^ 

For the additive Gaussian noise channel described by (16), 
each y(t.) is Gaussian distributed; therefore, this weighted 
sum of y(t^) samples must also be Gaussian. The mean and 
variance of G are easily obtained. Replacing x in (16) by 
and computing the expected value of G yields 
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E{G) - e[| J^x(tj-T„)[Ax(tj-T„) + o(tj)]j 


AP 


x‘ 


Also, the variance is easily found as follows. 


e{Cg-e(G)3^> 






Since E{n(t.)n(t )} = 0 for i ^ j , 

X J 


E{[G-E(G)]^} = K 

N 


I X (t -t) a 
i=l 


2 

n 


,2 r 

n 


N 


1 N 

- y 

N iii 


x2(t^-x)| 
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It is useful to define a normalized statistic with unity 
variance. Therefore, let 


and note 


g' = ^ G 


A 


E{G } 


«2 . 1 
0 , » 1 


It then follows that (G ) is a non-central chi-square 
random variable with non-centrality parameter [46] 


E{G 


NP^A‘ 


X - N(SNH). 


Note that (SNR) can be written as 


(SNR)= 




P (O ) 
X n 
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.2 


2 

P (a ) 
x' 


.2 

= G 

r N ■) "2. ■ 

W""' 

Also, from (22) and (23), 


Wow, fron (18) and (19) It is clear that S and (a^) 
represent the square of a weighted sample mean and the 
saiqple variance, respectively, of a Gaussian population. 
The weighted sample mean and the sample variance of a set 
of samples from a Gaussian population are independent [47]. 
Hence, the distribution for (SNR) is the ratio of two in- 
dependent random variables. Using the form in (23), the 
numerator has been shown to have a non-central chi-square 
density. The denominator, which is just a scale factor 
times the sample variance, is a central chi-square statis- 
tic with (N-1) degrees of freedom [48]. 

It is well known that the ratio of two central chi- 
square random variables leads to an F-distribution ; the 
ratio of a non-central chi-square random variable to a cen- 
tral chi-square random variable follows a non-central F- 
distribution [49]. If the chi-square statistics are 


(23) 


(24) 


normalized by their respective nu'^bers of degrees of 
freedom, a new random variable is defined using (23) by 
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.2 


F 


= (G )/»n 

r XT ^ 


f^l 

'^2 

(o2)/n 

L'^nJ 



(25) 


where 


m = number of degrees of freedom of the non-central 
chi-square numerator 

n * number of degrees of freedom of the central chi- 
square denominator. 

This ratio of normalized chi-square random variables has 
a non-central F density given [50], [51] by 

f (p’lm.n.X) * J e ^ V?’ |nH-2i . n) 

' i»0 ^ * 


with 


fp(F |k,n) 


k 

= (-) 

r(|)r(f) “ 


k-2 


IZ_1 


• • 

k ' 

1 + - F 
n 


k+n 

2 


= 0 


for 

f' >0 

for 

f' <0 


(26) 


where X = non-centrality parameter of the numerator of 
(25); it is given by (21). Finally, note that from (25) 


.2 




TN 1,"2. 
o 

n J 


(^) 

'm' 


F 
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so that using (23) 


F* = H (SNR) (27) 

m 

I ^ 

which provides the relationship between F and (SNR). 

Equations (26) and (27) determine the density func- 
tion for (SNR). For ease of computation, several approxi- 
mations to the non-central F cumulative distribution func- 
tion are available. One method involves approximating the 
non-central F-distribution with a central F-distribution , 
which is a widely tabulated function [52]. Another approx- 
imation given in the same reference makes use of the 
standard normal distribution function. This approxima- 
tion for the cumulative distribution function is given by 

F ,(f’ lm,n,X)=P(x ) 

F 


where 


r 0 ^ r 

mF 3 , 2 

Lra+X. 1 9nJ 

2(m+2X^ 

1 9(m+xy2j 

2 I m+2X \ 2 



and P( * ) is defined by 



(28) 


Substituting the values m*l, n=N-l, and X=N(SNR) from 
(21) into this expression yields an approximation for 
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the cumulative distribution function, F(SNR), in terms 
of the number of data points used, N, and the true SNR. 


F(SNR) =P< 


'(N-D(sk)l ^ r, _ 2 1 [ _ 2^-^lN(SNR) 1 

1+N(SNR) J [ 9(N-1)J [ 9( h-N(SNR))^J 


2+4N(SNR) 


9(1+N(SNR) )‘ 


9( 


2 f(N-l)i.SNRJ 
N-1) [l+N(SNR) 


0-1 


(29) 


Using (29) with P( • ) defined in (28), confidence intervals 
can be constructed. Ninety-nine percent confidence in- 
terval curves have been conq[)Uted for SNR from 0 to 20 dB 
asstaniG^ N*1000, N-5000, and N»50,000 sauries used for 

A 

the (SNR) estimate . Plots of the confidence interval 
curves for the white Gaussian noise channel are shown for 
SNR values from 0 to 10 dB in Figure 2 and for SNR values 
from 10 to 20 dB in Figure 3. For the values of N con- 
sidered, the widths of the confidence intervals vary less 
than 0.01 dB for SNR values from 20 dB to 40 dB. Hence 
the results at 20 dB may be extended to include most higher 
values of the SNR that might be needed. 

B. METHOD OF IMPLEMENTATION 

In order to form the desired estimates for the SNR and 
the normalized mean-square error, it is first necessary 
to find estimates for P^^, P^, and Perhaps the 

most straightforward method for obtaining these quantities 
is shown in Figure 4. The estimates for P^^ and P^ are 
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Figure 2. 


Ninety-nine Percent Conf dence Interval 
Curves. SNR = 0 to 10 dB 




Figure 3. Ninety-nine Percent Confidence Interval 
Curves, SNR = 10 to 20 dB 
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Figure 4. Direct Computation of P^, P^ , and 
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found by simply averaging the squared sample values. 

The cross-correlation estimate is found by delaying the 

samples from x(t), forming the products of samples 

x( c )y ( t ^ ) and averaging over the available .sampl(> 

pairs for any given delay i. Then Rjjy^'m'* found by 

computing R (t) for all values of t whert’ data ar<“ a\ a i 1 - 
xy 

able and picking Rv„('rm) to be the maximum valu(^ com- 

y 111 

puted. This provides the estimate for also. It is 

clear from Table I that the values for P^, and R^y^'m^ 

can be used to form all the desired estimates. 

For estimating P and P , this method is easily 

X y 

applied. Unfortunately, however, this direct approach 

for estimating R„„(t_) is extremely time consuming for 

xy *1* 

large sets of data. The entire summation of products 

must be performed for each value of t . A faster approach 

is to utilize the fast Fourier transform (FFT) algorithm 

to find R (t) by applying frequency domain techniques, 
xy 

and this is the procedure used here. 


The first step implemented in the computer routine is 
to find the approximate value of system delay u.sing the 
FFT. The necessary cross-correlation function, 
is the inverse Fourier transform of the cross- power spec- 
tral density, S (f). Thus 
xy 



(T) 


F‘Hs^^(f)] 

■xy 


= F~^[X*(f )Y(f )] 


( 30 ) 
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where the asterisk denotes the complex conjugate, and 

X(f) and Y(f) are the Fourier transforms of x(t) and 

y(t), respectively. Hence, it is possible to find the 

cross-correlation function between two sets of data 

x(t) and y(t) by finding their Fourier transforms X(f) 

and Y(f), conjugating X(f). forming the product 

X*(f)Y(f), and finally inverse Fourier transforming to 

obtain As explained earlier, the value of x for 

which R (t) is maximized is the value which the FFT in- 
xy 

dicates will minimize the mean-square error. Figure 4 
c«B be easily modified to utilize the FFT for finding 
Bjjy(Tjjj); the resulting configuration is shown in Figure 
5. 

It is convenient to be able to specify the size of 
the array to be Fourier transformed Independent of the 
total amount of data collected. Often it is helpful to 
utilize relatively short transform lengths. One reason 
is that errors due to roundoff increase as the size of 
the transform is increased [53j. Perhaps more important, 
the amount of memory required to be in the main storage 
of the computer is governed chiefly by the size of the 
transforms used. The approach chosen allows the FFT size 
to be specified as any value less than or equal to the 
total number of data values available. Of course, many 
algorithms for performing the FFT impose the additional 
restriction that the number of points transformed be a 
power of two [54]. 


Samples of x(t) 


Samples of y(t) 



P 


X 


R (V) 
xy m' 


P 
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Computation of P , P , and R„„(^) 
the FFT X y xy 


Figure 5. 


Utilizing 
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The X and Y arrays contain the sample values from 

x(t) and y(t), respectively. Before Fourier transforming 

the data, the X and Y arrays, which are of length K, 

are divided into smaller blocks of length N. This is 

illustrated in Figure 6. Because of the assumed wide- 

sense stationarity of x(t) and y(t), it is possible to 

estimate R (t) by selecting a short block of data from 
xy 

the X array and correlating it with the appropriate 
sauries from the Y array for all values of t. The tech- 
nique used to find Rjjy(f) is to transform blocks of data 
of length n one mt a time and form Rjjy(t) as a nui^r of 
separate segments. The block of | samples taken from the 
beginning of the X array is padded with zeroes to form a 
block R sauries long. This serves as a reference block. 
Then the first block of N sain>les is taken from the Y 
array. The reference block and the first measurement data 
block are Fourier transformed, and a segment of R^y^^^ 
found using the method described. Of the N values which 
result, the first | values are valid values of 

The last — values are extraneous data and are discarded. 

2 

The last values fail to be valid because of the periodic- 
ity assumed in using the FFT algorithm. This cyclic 
property of the FFT is also the reason for padding the 
reference block with zeroes. 

In order to find the next ^ values of R^y^^^' ^ 
block of N samples is taken from y(t) which overlaps 
exactly one half the first block. This is also shown in 
Figure 6. These samples are processed just as lor the 
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first block. The process Is continued ss long as there 
are data points In the Y array to transforr,. i„ this way 
Ray(^) is found for all delay values, except tor the very 
large values where there .s Insutticlent data to pertor™ 
the last PFT. The array is then scanned to find 

the delay value for which it is a maximum. Pinding 
«xy<f) in this way allows the length of the transformed 
arrays to be Independent of the total amount of data to 
be analyzed, except that total data length serves as an 
upper bound on transform length. 

once a delay value has been found using the PFT a 
direct tl« do«i„ approach is used to compute seve^l 
TOloes of h^^(t) in the vicinity of the delay Indicated 
by the PFT. There are two major reasons for doing this 
0~ is that the PFT IS affected by cumulative roundoff 
t^rors to a greater extent than is the direct (time 
domain, method. Another is the fact that each value of 
the correlation function computed with the PFT is obtained 
using only 5 data values. The time domain method utilizes 
all the data available from the current delay value to 
the end of the data stream. Additionally, the time domain 
calculation includes normalization by the power in the 
samples being used at each Iteration of the delay. This 
corrects for small changes in the power contained in the 
array as single samples are added and deleted from the 
Rxy(t) calculation. After the delay is found, the other 
quantities follow easily. Provision is also made to use 
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a Ham&ing data window, if desired, to smooth the final 
estimates for Rj^y(t), P^. and P^. 

A block diagram of the complete routine is shown in 
Figure 7. Input data for the routine consists of the X 
array of reference data, the Y array of measurement data, 
and various control parameters which indicate such things 
as whether to window the data, what the FFT block length 
is to be, etc. As shown in the diagram, the reference 
block is selected from the X array and loaded into a work 
array XT for transforming. A block of measurement data is 
selected from the Y array and loaded into the work array 

TT for tr^sldfalBg. Botft XT and YT are traBsformed by 

' ‘ # 

the FFT, aad the product X (f>Y<f > is inverse transforraed 

to yield R (t). This process is repeated for each block 
xy 

of measurement data, and the delay value where the FFT 
indicates a peak in is found. Then iterations are 

performed for a few values of t around this peak, and the 
direct time domain method with double precision arithmetic 
is used to compute 


P 


»^P~P~ 

X V 


for each iteration, where 


UU) 




1 

K-t 


K-t 

I X(i)Y(i+T) 
i»l 


( 32 ) 


and 
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7. Block Diagram of SNRMSE 


Figure 
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K--t 


= -1- "y' x^(i) 

1/ -r ^ 


( 33 ) 


K-t 


i=l 


and 


K 


= — I 

ir 'T ^ 


k-t 


Y^(i ) 


( 34 ) 


i=i + l 


K i«? the total number of data 
In the above expressions, K i 

swnple index In the array. The nun,ber of iterations per- 
formed about the initial estimate lor t Is specified hy 
an input parameter, ItHOAII. The peak of this nommllaed 
cross-correlation defines The value of Py is computed 

directly as the average of the squares of the samples in 
Y from (t„tl) through the end of the data. The value of 
Is similarly computed using the first (K-t„) samples 
In X. The Hamming window is applied to the data, if 
specified, in doing the final computations. These same 

sets of data samples arc used to compute the cross- 

V R ^ It is important in calculatinR 

correlation peak, Rj^y(T^)* 

P and P to use only those samples which were used in 

finding estimate for system gam should 

. . _ Q ^ as it is in ntcasuring 

be the same in measuring m 

p In computing the SNR, for example, using 

y' 


- "xy'^m’ 


SNR = 
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the quantities (P^Py) and are Identical in the 

first three, four, or even five digits for a high SNR. 
Thus it is critical to have identical estimates of gain 
implicit in these two terms. In order to make the mean 
square error measurements independent of the effects 
of system gain and power in the reference data, a 
normalized mean-square error, BURN, is computed using 


ERRN= — B 




( 35 ) 


Double precision arithmetic is used for computing all 
final output values. A listing of the computer sub- 
routine SNBMSE is given in Appendix A. 

There are several things which must be considered 
in using the computer routine for SNR estimation if 
accurate results are to be obtained. Two obvious con- 
siderations concern the total number of data points used. 
Clearly, more data points will usually provide more 
accurate results. Also, the data stream must be long 
enough to make certain that the system delay value falls 
well within the total length of the Y data array being 
analyzed. Choosing a size for the FFT block involves 
more subtle considerations. Errors in the FFT algorithm 
increase as the block length of the transform increases, 
which would indicate smaller block lengths should be 
best However, the use of small block lengths forces 
select the preliminary delay value (the 


the routine to 
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value indicated by FFT) based upon short data segments. 

Hence, the routine is more likely to miss t completely 

m 

if the block length is small. One key assumption in 
the development of the routine is that the data is 
stationary from block to block. If the block length 
is chosen to be very short, this assumption may fail to 
be valid. The routine tests the mean and variance of 
each block to be transformed and generates a warning 
statement if the data does not appear approximately 
stationary. A brief description of these stationarity 
tests Is given In Appendix A. These considerations In- 
dicate that very short FFT block lengths are also un- 
desirable. ScMse typical values of block length which 
have been used with success are 2048, 4096, and 8192. 

All these values of block length seem to give excellent 
results. It does not appear that the largest block 
lengths used significantly degrade accuracy through round- 
off errors. Block lengths less than 2048, however, have 
occasionally been found to miss the initial delay 
estimate by large amounts. Results of tests using various 
block lengths are given in the section on applications to 
analog systems. 

It has also been seen that the window function pro- 
vides no noticeable improvement in the accuracy of the 
routine with most of the systems tested. For this rea- 
son, windowing was seldom used in the final testing of 
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the lowpass SNR estimator. It .will be seen that this 
is not the case, however, for the bandpass SNR estimator. 
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V. AN SNR ESTIMATOR FOR 
COMPUTER SIMULATIONS OF BANDPASS SYSTEMS 

A. THEORETICAL DEVELOPMENT 

The problem of estimating the SNR in simulations of 
bandpass systems is complicated by the fact that simple 
time-shifting and amplitude scaling of the reference will 
not yield a replica of the signal component of the mea- 
surement data, as it does for lowpass systems. Estimates 
for both the grotq> and phase delays must be simultaneously 
determined. This probl^ is solved by confuting all the 
required quantities in the frequency domain. Subroutine 
SNRBPS is a computer routine which was developed to 
estimate the SNR and other parameters in simulations of 
wide-sense stationary bandpass systems. 

As in the lowpass case, a reference signal x(t) and 
a measurement signal y(t) are required. The estimate for 
the signal component of y(t), denoted by z(t), is defined 
to be identical to the reference signal except for some 
unknown amplitude scaling and unknown group and phase 
delays. To form a minimum mean-square estimate of the 
signal, it is necessary to define z(t) such that the error 

E = [y(t)-z(t)]^ (36) 

is minimized, where the overbar denotes a time average. 
This implies that the power on the given interval in 
the difference signal [y(t)-z(t)] is to be minimized. 
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As a direct result of Parseval's theorem t55], it Is 
clear that this Is accomplished If power computed from 
the spectrum of Cy(t)-z(t)] on the Interval Is minimized 
The FFT may be used to compute the desired sped 
and the error can be expressed by 


N 

= I 

i=l 


Y(f^) 


- Z(fi) 


vrhere Y(f^) and Z{f^) denote the ith samples from the 
Fourier transforms of y(t) and z(t), and N is used 
denote the nun^er of points transformed. For convenience, 
it will be assumed that N is even; if not, it can easily 
be made so by deleting one of the original data points. 

No significant accuracy is lost since only rather large 

data sets (N>500) are likely to be used here. If a 

wss-wn criterion is automatically 

binary radix FFT is used, this crixerxw 

satisfied already. The symmetry Involved In the FFT 
algorithm makes It clear that the summation over the last 
I complex frequency samples must be identical to that over 
the first I samples, except lor errors within the al- 
gorithm itself such as roundoff or truncation. This 
analysis will neglect such errors. In the actual Imple- 
mentation It will be explained later that all samples 
are utilized in an effort to minimize roundoff errors. 

For now. however. It will be expedient to consider mini- 
mizing the error computed from only | samples. Hence the 
error expression to be minimized can be 
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N 2 

1 2 

i E = I Y(f . ) - Z(f . ) 

^ 1=1 ^ ^ 

which represents one half tht' total error power. 

Consider some al^obraic manipulations on the 
error expression as follows: 

N 

1 2 , 
i E = I [Y(f.) - Z(f.)][Y(f.) - Z(f.)] 

-s i^l 1 11 1 

N 

2 * * ♦ 

= [Y(f^)Y (f^) + Z(f^)Z (r^)-Z (f^)Y(fp 

- Z(f.)Y*(f.)l 

N 

2 

= I { |Y(f )|^|Z{f )!^-[Z*(f . )Y(f , ) + (Z*(f . )Y(f . ))*'} 

i=l ^ 1 11 11 

N N N 

2 2 2 

= I |Y(f )|^+ I jZ(f )|^- y {2Re:Z*(f )Y{ f ) } 
i=l ^ i=l ^ i=l ^ ^ 

where Re[«] denotes the real part of a oomple.x quantity. 
Thus, the error expression to be minimized can be written 

N <? N r N 

2 I I 2 2 

I E = I lY(f )| -I- I |z(f )P - 2Re I Z*(f )Y{f ) 

i=l i=ll ^ ' Li=l 
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Since Z(f^) is not directly available, it is neces- 
sary to substitute 9 representation of Z(f.) in terms of 
the reference X(f^). From the definition of z(t) given 
earlier, it is clear that Z(fj^) should ditfer t rom 
X(f.) only by the frequency domain equivalents of group 
delay, phase delay, and gain. It follows from the defi- 
nitions of group and phase delay [56] that, for the 
assumed linear phase characteristic, phase delay of the 
signal estimate z(t) can be adjusted by adding a constant 
fdtase wgle to all frequency coa^nents. Also, the 
grotqi delay of the estimate for the signal can be ad- 
justed by adding a phase angle to each ctw^jonent which 
is directly proportional to the frequency of that com- 
ponent. Gain is adjusted by a scale factor. Just as In 
the time domain. Therefore, the expression for Z (f.) 
for positive frequencies (the first | samples from the 
FFT) can b© written 

J(af.+b) 

Z (fj) = AX (f^)e (37) 

where A, a, and b are real constants which can be chosen 
to minimize error. For negative frequencies the sign of 
b is reversed. Adjusting the parameter b is equivalent 
to adjusting the pha^e delay estimate in z(t). Varying 
the parameter a provides for adjustment of the group 
delay estimate. 
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Substituting (37) into the error expression to be 
minimized yields 


N N 

2 I 2 2 -j(af^+b) 

I Y(f^ I AX(f.)e 

i*l| i=l 


r N 

2 * j(af.+b) 

I Y(f.)AX (f.)e 
i«l ^ -» 

w 


- 2Re 


If E is used to denote the minimized error expression. 


2 g 


5 1 -J(atj+b)r 2, 

E« - MIR S ! AX(f.)e + I|Y(fi)| 


A,a,bji»lf 


r 

2 * j(af.+b) 

- 2Re I Y(f.)AX (f^)e 


MIN <A^ I lX(f.)r-^ I 17(1^)! 
,a,b| i*l i=l 


2ARe I Y(f.)X (f^e 
i=l 


Jaf, 


i 


(38) 


It should be remembered that E« equals one half the total 


mean-square error. 
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Clearly the first two terms in (38) are positive 
and independent of a and b. Thus, assuming that A is 
a positive scale factor, the optimal choices for a and 
b must be such that 


Re 


r N 

^ lb 

I Y(fpX (ip e 

i = l 


is maximized. Therefore, consider the expression 

N 




✓ 

/ 

Re 



\ 


{ ! *(fi>** 
\i-i ‘ 


(fi>e 


Jaf 



( 39 ) 


Note that once the summation has been performed, the 
resulting sum will be some complex number, say C(a). 
Then (39) becomes 


MAX{Re[e*^^C(a)]l . 
a.b 

If the angle associated with C(a) is denoted by 0, then 
C(a) = lC(a)|e*^®. Multiplication of C(a) by yields 

e-^^C(a) = e-^^^'*’®^lC(a) I . 

The vector represertation for C(a) and e'^^*C(a) are 
illustrated in Figure 8. 

The angle of [e^*^C(a)3 may be adjusted to any 
desired value by varying b, and the magnitude, le'^^C(a)l 
= lC(a)l , is independent of this adjustment. It is 
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clear from the diagram that the real part of [e^*'C(a)] 
is maximized when b » -e so that [eJ*’C(a)] Is real and 
positive. For this choice of b, 

Re[e-^'^C(a)] = lC(a)l . 

Therefore, the maximization over b always leads to 


MAX{Re[e^^C(a)]} = MAXlC(a)|. 
a,b a 


(40) 


Althougb this developnent is valid, it should probably 

be mentioned that the derivation of equation (40) is 

simplified somewhat by the decision which was made to 

utilize only | samples in the theoretical development. 

If all the samples are used, the summation over N samples 

, - . N 

must be broken into two sums over the first and last ^ 
samples, and discussions about the symmetry of the FFT 
are required. After some development, however, a 
criterion for choosing b can be obtained that is the 

same as (40). 

Using the result given in (40). expression (39) may 
be written more simply as 


MAX 

a 


I Y(f . )X*(f -)e' 
i=l ^ 


jaf 


(41) 


With a = 2 ttt the preceding expression may be rewritten as 
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N 



or 


N 

I 2 

MAX I I 

T i=l 


S (f . 
xy 


j2Tlf T 

)e ^ 


(42) 


where is the ith sample of the cross-spectral 

density [X*(l)va)], Let t„ denote the choice ol t which 

n.«lmlzes (42). The the value of Is the group delay 
estioKtew 


Once has been found. It is useful to define 



(43) 


Then , 


(38) and (43) can be used to obtain 



Denoting 


(44) 


N 

1 2 ^ 

2 power in X = P = ^ |X(f.)|2 

1=1 ^ 

and 

N 

1 2 „ 

2 power in Y = f |Y(f.)|2 

i=l ^ 


provides the error 


expression 



67 


Eg = MIN[A Pjj+Py-2AR2] 
A 


To find the value of A that minimizes Eg, the deriv- 
ative with respect to A is computed and set equal to zero, 
This yields 


dE„ , n 

asr - as [aV’y-2*"2>® 


Using to denote the optimal choice for A, 




Thus, 


^2 

Ajjj = ^ = System Gain. 


(45) 


Substitution of (45) into the error expression 
yields 


^2 = 


A P^<-P„-2A R„ 
m X Y m2 



or 


^2 = 


(46) 
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The signal power in the measurement data is clearly 


S = A^(2P„). 
m A 


(47) 


From (46) and (47), and remembering that E2 represents 
one half the total error power, the signal-to-noise 
ratio is given by 


SNR = 


S 

2E, 


II 


^m^X 


P -A^P 


or in terms of Rg, Pj, and Py, 


where 


SNR 


W 


r: 


N 

2 


P„ = 


I |X(f^)l 

i=l ^ 


P 


Y 


N 

2 


I lY(fi)l 


2 


and 


N 

2 * 

I X (f.)Y(f ) 
i=l ^ 


e 


(48) 
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Estimates for various system parameters that are avail- 
able in terms of these quantities are tabulated for 
reference in Table II. 

As in the lowpass case, it can be shown that the 
bandpass SNR estimator provides the ML estimate lor the 
additive white Gaussian noise channel. In this case a 
frequency domain approach is used. It is necessary to 
determine the joint density of the frequency domain samples, 
Y(f^), obtained by the FFT, in terms of the gain A, the 
group delay parameter a, and the phase delay parameter b. 

First note that for the additive white Gaussian noise 
cbannel y(t) ■ z(t) n(t). TOerefore, the FFT yields 

Y(fj) « Z(f^) + N(fj) (49) 

where Y(f^) denotes the ith frequency component of the 
waveform under test, y(t), and Z(f^) and N(f^) represent 
estimates for the signal and noise components, respectively. 
The discrete transformation of the noise may be written 

, N-1 -j2TTf.k/N 

N(fi> = 4 I n(t ) e 
^ ^ k=0 ^ 

or 

1 k 

N(fi> = N I n(tj^)cos(2irf 

0 

‘ ^ ® Jo S> 



TABLE II 

PARAMETER ESTIMATES FOR BANDPASS SYSTEMS 
IN TERMS OF Pjj, Py . ^2 
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which Is 

N(fi> = Njj(f.) - 3 

where 

H Jo n(tj^)cos(2iTf . |) 

and 

^ J^n(tj^)sin(2itf . |). 

The parameter N represents the total number of points 
heli^ trassforBed. 

Since and linear coBd>inations of 

Gaussian random variables, they must themselves be Gaussian 
distributed. Also, they represent the direct and quad- 
rature components of the 1th spectral component and are 
uncorrelated [57]. Thus for any given i, Nj^(f.) and 

are independent Gaussian random variables. Further, 
since spectral components of different (non-overlapping) 
frequencies are uncorrelated, it follows for any ergodic 
process that Nj^(f . ) and Nj^Cf^) must be uncorrelated for 
any i j where the FFT gives a valid representation of 
the spectral components in the frequency domain [58]. 

The same argument holds for Nj(fj^) and Nj(fj). The FFT 
provides a valid representation of spectral components 
for frequencies up to one half the sampling frequency , 
if properly applied. This corresponds to the first 
half of all the computed frequency domain samples. 
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signals, which ar© th© only on©s consid©rod 
N 

here, the latter g complex samples produced by the FFT 
are the complex conjugates of the first ^ samples [59]. 
This deterministic symmetry makes it clear that no addi- 
tional information is gained in using the second half of 

the samples. Therefore, the analysis continues based 
N 

upon only ^ complex samples. 

Uncorrelated Gaussian random variables are inde- 
pendent. Therefore, the first half of all the noise 
samples computed by the FFT. Nj^(f^) and Nj(f.), are 
independent. 

It is useful to extend the R and I subscript notation 
to include Y(f^) and Z(f^) so that 

(real components) 

and 

Yi(f^) = Zj(f^) + Nj(f^) (imaginary components). 

Then, the joint density of the first | real sample 

values and ^ imaginary sample values can be written as 

N 

2 

P_(Y|A,a,b) = Pi(Yjj(f^)|A,a,b)p.(Yj(f . )|A,a,b) 

where p . ( Yj^( f . ) | A, a ,b) and p j (Yj ( f ^ ) | A , a , b ) are the 
Gaussian densities for the ith real and imaginary sample 
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values, respectively. Since N(fj) Is zero mean, the 
means of V^tfj) and ViCtj) are given by Zp(fj) and Zjdj). 
Under the white noise assumption, the noise variance will 
be the same for each N (fj) and and will be denoted 


by 0 ^. Admittedly, this is not generally the case, since 
filtering is usually present which scales noise dif- 
ferently at different frequencies. However, an exact 
analysis would be system dependent and a general solution 

does not appear possible. 


From (37) it is seen that 


» A[Xjj(f^)cos(af^-^b)-»-Xj(f^)sin(af j^+b)] 

and 

Zi(fi> » A[Xj(f j^)cos(af^+b)-Xjj(f j^)sltt(af^+b)]. 

Using these results and assumptions, the joint density 


becomes 

N 

H 1 

P (YjA.a.b) = n 2 

Y i=l 



74 


This density, which is the likelihood function ior A, 
a, and b can be rewritten after some algebraic manipula- 
tion as 


P_(Y|A,a,b) = 
Y 



( 52 ) 


The log-likelihood function is 

2 * j(af +b)L 

2ARe y X (f )Y(f . )e 

i=l i " J 


N 


2o 


N 


N 

2 

- I 

i=l 


Y^(f . )+Y^(fp 


I 


i = l 






-|ln(2TTc^j) 




(53) 


By definition, ML estimates for the parameters a and b 

are those values which maximize this function. Assuming 

2 . 

that A is any positive scale factor and is any positive 
value, the optimal choice for a and b must be such that 
the expression 
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Re 


r N 



Li = l 


j(af +b) 
(f.)Y(f.)e ^ 


( 54 ) 


is maximized. Note that this is true regardless of the 

2 

values of A and which will be variables in determining 
signal and noise estimates later. This leads immediately 
to the criterion expressed in equation (39), and thus it 
is seen that the estimates formed for parameters a and 
b are ML estimates. Parameters a and b are directly pro- 
portional to the group and phase delay estimates, res- 
pectively. Thus the group and phase delay estimates 
obtained with the minimum mean-square error criterion 
are UL estimates for the additive white Gaussian noise 
case. 

Once a and b have been chosen, expression (54) becomes 
a constant. Denoting this constant by , the log-like- 
lihood function becomes 


N M 



N 2 

- I ln(2Txap. 

At this point it is convenient to drop the subscripted 
form of X and Y. Rewriting in terms of the original 
complex form yields the log likelihood function 
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^ i = l 
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i = l 


X(f.) 


I ln(2 
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From (43) and its earlier equivalent form. (39), and the 
definition of , expression (54) maximized over a and 
b. it is seen that R, = R«. Usins this result and the 
definitions for and Py given following (44), the 
log-likelihood function can be written 


- I ln(2ir0j). 


Just as in the lowpass case, the UL estimate for the SNR 

is formed as the ratio of the ML estimates for the signal 

power and for the noise power. As a convenience in using 

the notation developed for the bandpass case, the estimates 

will actually be formed for one half the signal power and 

one half the noise power; obviously the ratio is un- 

2 

changed. One half the signal power is given by A P^. 

2 

Defining ^=A P^^, remembering that A is assumed to be 
positive, and rewriting the log-likelihood function 
yields 


2o 


N 


2 - P„ - 

2 ^ Y 2 

X 


N 1 /o 2, 

- 2 ln(2i!Ojj) 
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One half the noise power is given by (N o^) so that the 
log-likelihood function can be written 


N 


2(Na^) 


2R 




^ - Py ^2 


-i- 


2ti , M ,2 V 

tr*'* “n* - 


(55) 


The ML estimates for S 2 and (N Oj^) can be found by 
setting the partial derivatives with respect to and 
(N 0 ^) equal to zero and solving. For the signal power, 
the partial derivative of (55) yields 


N 

2(N 0^) 




1 


= 0 


or 



(56) 


For the noise power, the partial derivative of (55) 
yields 




2R„ pc— 

K - Py - 


(N 5^)^ 


N 

2 


(N ap 


= 0 . 


Substituting from (56) and simplifying 
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''2 

(N of^) 


+ 


1 = 
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or 


(N 




(57) 


As explained in the development for the lowpass 
the ML estimate for the SNR is simply the ratio 
ML estimates for the signal power and the noise 


eas<' . 
of t h»' 
powt'r . 


Thus 


SNR 



2(N aj) 


(N op 



or 



(58) 


The estimate given by (58) is precisely that of (48) 
Thus for the white additive Gaussian noise case, the 
estimator developed is the ML estimator. 


Unfortunately the development of confidence interva 

for the bandpass case does not appear to be mathematical 

o 

feasible. The computation of the quantity involves 

the sum of the squares of two Gaussian random variables 
which are not easily shown to be independent. This may- 
be seen by considering 


Is 

ly 


N 
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Rp = IX (f )Y(f )e 
i = l 


j2irf .T 
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i=l 


I [XR(fi)-JXiCfj)][YR(fi)*JVi(fi)][cos(2i.t,T^) 

+j sin(2Ttl-.T|^) j 


r N 
2 


j;^R(fi)[Xj^(f . )cos( 2 TTf .T^)+Xi(fi)sin(2iTf .T^) j 

- Y-(f^)[XR(l.)sin(2wf .T^)-Xj(f^)cos(2Tif^T^)]j 


f E 
2 


I YR(l^)[XR(fj^)sin(2wl .Tjjj) - Xj(f ^)cos(27rf 

i=lL 


+ Yj(f . )[XR(f .)cos(2TTf . '^)+Xj(f . )sin(2TTf -Tj^)] 


An inspection of the various terms within braces indicates 
that it may be possible to show the individual terms 
within the first squared quantity to be independent of 
the individual terms within the second squared quantity. 
However, this is insufficient to guarantee independence 
of the complete quantities. If independence were assumed 
in order to continue the analysis, the distribution of 
rJ would be non-central chi-square. It does not appear 
that R? is the square of a weighted sample mean of a 
Gaussian population. This means that the denominator 
of (58) cannot easily be shown to be the sample variance 
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of a Gaussian population. Hence the distribution for the 
denominator of (58) does not follow easily, as was the 
case with the lowpass estimator. Finally, it is not 
clear that the numerator and denominator in ecjuation 
(o8) are independent. Thus, determination of the dis- 
tribution function for the bandpass SNR estimator appears 
to be a very formidable problem. It may be FK)ssible to 
form an approximation to the distribution function, how- 
ever. The estimate has been shown to be the ML estimate 
for a white additive Gaussian noise channel. It can be 
shown that under reasonably general conditions, the ML 
estimate Is assnnptotically Gaussian as the number of 
samples used approaches infinity [60]. Thus, one approach 
might be to develop a Gaussian approximation, especially 
since the cases of interest generally involve at least 
1000 samples. The difficulty with obtaining the approxi- 
mation arises in attempting to determine the variance 

>v 

of (SNR), which is necessary for forming the Gaussian 
density. This does not appear to be a simple matter. 
Therefore, confidence intervals were not derived in this 
work for the bandpass estimator. 

Confidence intervals would be convenient in order to 
precisely define the reliability of the estimator, at 
least for the additive white Gaussian noise channel. A 
deeper statistical study of the reliability of the estima- 
tors IS an area in which further work could certainly be 
iiowevei , a good case for assuming adequate 


done . 
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reliability of the bandpass estimation routine can be 
made based upon empirical data. Results of tests of this 
estimator appear later in the section on applications to 
analog systems and throughout the section on applications 
to digital systems. Excellent reliability is evident in 
these results. 

B. METHOD OF IMPLEMENTATION 

In order to estimate the SNR and other parameters 
of interest in simulations of bandpass systems, it is 
convenient to first find estimates for Py, and Rg. 

As can be seen from Table II, a variety of system param- 
eters, including the SNR, can be computed in terms of 
these quantities. 

From the definitions given for P^^, Py. and Rg follow- 
ing (48), it is evident that the computation of the re- 
quired quantities is done after the FFT has been applied 
to the reference data, x(t), and the measurement data, 
y(t). A conceptually direct method for computing Py, 

and R„ is shown in Figure 9. However, as in the lowpass 
case, the most straightforward method of implementation 
is not practical because it requires that the complete 
summation over all complex samples be repeated for every 
possible value of t in order to find and thus determine 
R . This operation requires a great deal more computer 
time than is practical using the number of samples 
generated in a typical simulation. In order to make the 
SNR estimation routine practical, some method must be 
found for rapidly determing 
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Consider the quantity to be maximized in locating 
r^. From (42) and (43) 


R„ = 


MAX 

t 


N 

o 

! 

i = l 


S (f.)t 

XV 1 


.l2nf . I 


where S ,(f.) = [X (f.)Y(f.)l and i is the valiu' of i 
X y 1 i i m 

for which maximization is achieved. Consider the complex 
quantity generated by the summation and define 


N 

2 J2irf T 

"a*’’ ■ Ji ® 


Also, make the definition 


S.(f . ) 

A' 1 ^ 


S (f . ) 

XV 1 


- ■ N 

for lli^ ^ 


= 0 


for (^ + 1 )l i i:N 


This yields 


N j2Tif.T 


Let T be chosen to be some value of time eciual to 
an integral multiple of the time interval between samples 


so that 
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T = kT, 


0<k<N-l 


where T is the time between samples. V/ith this 

this restric— 

ion R^(t) is, by definition [611 thc> • 

_ inverse discrete 

ou„er transfer, of s^(f,. frec,ua„cy domain samples 

or , Wit, eanoe Viewed as samples fro. t.e 

- ive fredueney ran.e Ceal. TPns, 

the inverse Fourier transform of s (f) „,th rh 
a . with the nega- 

frequency components set to zero R-m k 
o ff. . zero. Remembering that 

xy(n IS the Fourier transform of R (x) it « 

D / w IS seen 

that R^(t) Is the function whioh », 

to K fir * ld«tlcal 

*y(t) for positive freouencies and equal to zero for 

negative frequencies. It follows incedlately from the 

definition Of analytic signals C63J that R,(t, .3 pro- 

portional to the analytic sie-nai 

^ytic signal corresponding to the 


real -Signal" R^^Cx). Then 


^2 MAX I R^ ( T ) 


and IS the value of x for which i-hi:a 

hich the magnitude of the 

analytic signal corresponding to R (,) ^ 

, xy ' ^ *^aximized. 

IS possible to generatf* thic 

generate this analytic signal quite 

rapidly with the FFT Thio 

IS approach provides the neces- 
sary technique for estimating gnlehly. 

It is possible to Obtain the desired analytic 
signal directly from x<t) and y(t, using the FFT This 
can be accomplisbed by forming r 2 X*(f,y(f,j 


85 


the last I sample values (corresponding to negative fre- 
quency samples) to zero before inverse transforming. 

An alternative approach is to compute R^y(i), whi«.h is 
the real component ot t ht' analytic signal , just as in 
the lowpass estimator and t<. find the imaginary component 
of the analytic signal by computing the Hilbert trans- 
form of Computation of the imaginary component 

may be done with a discrete Hilbert transform algorithm 
or with the FFT, using the simple frequency domain equiv- 
alent of the Hilbert transform [64], [65]. Since the 
software to compute and store been written 

and thoroughly tested in the lowpass SNR estimator, the 
latter approach for computing the analytic signal was 
used, and the Hilbert transform was performed with the 


FFT. 

Instead of computing the Hilbert transform and de- 
termining the analytic signal for the entire range ot 
R .(t), a preliminary scan of H^y ^ ^ performed, and 

the peak value of O is found. It is assumed that 

although the maximum magnitude of the analytic signal 
does not necessarily correspond to the peak of ^ ’ 

it is highly probable, using the slowly varying envelope 
concept, that the two maxima occur in the same neighbor- 
hooi. Also, it was experimentally determined that 
approximately twenty well behaved artificial sample 
values added at each end of the FFT transform block 
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helps to avoid any problems with errors resultinp from 
the discontinuities which may exist there. 

The twenty artificial values at the beginning of 
the transform block were obtained by taking the first 
twenty true sample values in reverse order . For example , 
suppose the artificial samples are denoted by YT(i), 

i*l 20, and the true data values follow, beginning 

with YT(21), YT(22), and so on. The values of the 
artificial samples are assigned using YT(21— i) — YT(20+i), 
foy i * 1 20. A similar procedure is utilized 

to develop tlie artiftcia saaq;»le values at the eud of ttie 
traBsfomn block uaiag the last twenty true data values. 
Experimental testing of this scheme showed that good re- 
sults were obtained for a variety of signals, so the 
method was adopted for use in SNRBPS , 

Therefore, in the implementation used to determine 
group delay, the assumption is made that the maximum 
magnitude of the analytic signal and the peak value of 
R (x) are separated by a delay value corresponding to 
no more than + - 20) samples, where N is the number 

of samples chosen for the FFT block size. Under this 
assumption only one block of N samples needs to be 
processed by the Hilbert transformation procedure, since 
the analytic signal is required only on the interval 
extending + (^ - 20) samples from the delay value where 
the peak of R^„(t) is located. 

Ay 
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This assumption saves computation time and avoids 

the problems associated with end effects that result in 

attempting to Hilbert transform many st'parale blocks of 

data to produce a longer transformed array. Ii the signals 

involved are not narrowband, it is possible that the p<*ak 

of R (t) could be sufficiently removed from the desired 
xy ' ^ 

N 

delay value that the assumption of a + (^ - 20 ) sample 
range could result in an error. This never occurred in 
any of the applications of the bandpass SNR estimation 
routine, however, and the assumption appears to be justi- 
fied by tt» resulting simplifications and reduction in 
computer time requirements. The possibility of error due 
to this assumption could, of course, be eliminated by 
Hilbert transfoming the entire array or by 

computing the analytic signal for the entire R^^Ct) array 
directly as mentioned earlier, if it were felt to be 
necessary . 

As with the lowpass SNR estimation routine, the 

elaborate processing performed to initially estimate 

delay serves only to determine a starting point for the 

final search. The successive Fourier transformations 

used to obtain R (x), and to implement the Hilbert trans- 

xy 

form, reduce the accuracy of the initial R2 estimate. 

Also, as in the lowpass case, the original R.,„(x) compu- 

Ay 

N 

tation for any given i value is based upon only 2 samples 
from the x(t) and y(t) arrays. As a result, the initial 
estimate obtained for R2 is unsuitable for forming an 
accurate estimate of the SNR. 
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A number of iterations are performed using a more 

direct and accurate computation to form trial values foi 

R„, denoted r 1 , in the vicinity of the initial delay 

estimate. The number of points to be checked in the 

final iterations is spec’itied by an input parumt’ter to 

the routine. Each of thest> computations is performed by 

first shifting the x(t) data to account lor tht' dt'lay 

vajue being considered and thus the exponential factor 

in the defining sum for 82* The delayed data are denoted 

bv Then the FFT is applied to data blocks from 

' d' ' ^2TTf -T 

mMmErnmmz ^ to , form. f j > * X( f i ^ 

and t{f respectively. Finally, complex conjugates 
are con^Juted, multiplication of the sample pairs is 
accomplished, and the suranation is performed. The raagni- 

S' 1.1 i ■•4'' ji'- it •i' "'•V ‘ V "• :'>■ "••■•! i". !v .-r ^ ^ ^ ^ • 

tude of the resulting sum is computed at each iteration 
and is normalized by /P^Py , where P^ and P^. are the powers 
computed from the same frequency domain samples that 
were used in computing R.^ - The peak of this normalized 
cross-spectral power computation defines and for 

use in the final SNR estimate. 

The reason for shifting the xCt) data to acctnmt 
for the factor before transforming instead of in- 

cluding such a factor in the frequency domain calculations 
for R^ is to avoid problems which occur because of the 
periodicity inherent in the FFT. Performing the frequt^ncy 
domain equivalent of time-shifting with a finite block 
length yields errors equivalent in the time domain to 
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overlapping the beginning of the block of data from one 
array with the end of a data block from another array in 
a cyclic fashion. This problem can be overcome by padding 
one array with zeroes before transforming as was done in 
initially omputing Rj^y(T) in both SNRMSE and SNRBPS . 

This still leaves some problems with high frequency com- 
ponents generated by the discontinuities which exist after 
abruptly setting the waveform to zero at some poxnt. Un- 
like the initial Rjjy(T) estimation, the result of the 
final Rg computation is used directly in the estimation 
of the SNR. Therefore, maximuiB possible accuracy is neces- 
sary. Time-shifting the array before transforming avoids 
most of the problems. Although slightly nrare complicated 
in concept, this approach is simpler in implementation 
and more accurate. 

Actually, the basic quantities which are computed 
in the final implementation are (2P„), (2P^), and (2R„). 

AY 2 

These computations result from the fact that P , P , 
and R 2 have been defined as stimmations over only the 
first half of the complex frequency domain samples com- 
puted by each application of the FFT. The symmetry which 
theoretically occurs when the FFT is applied to real data 
indicates that including the second half of the complex 
samples would yield identical results to those obtained 
using only the first half. In practice, this is not 
quite true because of roundoff errors which occur within 
the algorithm due to finite register length. It seems 
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unlikely that roundoff errors which appear in the 
second half of the samples computed by the FFT would be 
identical to th('>se in the first half. Also, it is logi- 
cal to assume that the mean-square value of the error 
introduced by roundoff is roughly the same for samples 
in both halves of the computed array. If these assump- 
tions are true, the effects of roundoff error should 
be reduced by utilizing all the complex frequency domain 
sanples, instead of oely half of them. The additional 
samples are readily available and the additional time 
reQuired to utilize them is not a significant factor in 
the overall computer time requirements of the estimation 
routine. Therefore, in an effort to minimize the effects 
of roundoff errors, all samples computed by the FFT are 
utilized. Summations over all the samples produce the 

quantities (2Pj^), (2Py). and (2R2)» 

In performing the final computations, a data window 
is applied before Fourier transforming the data to smooth 
the frequency domain estimates, if specified using an 
input parameter to the routine. A variety of data win- 
dows were tested for this purpose. It was found that 
the application of data windows often improves the SNR 
estimates in cases where the SNR is relatively high. The 
Hanning window appeared to yield the best results and 
was chosen for use in the SNR estimation routine. Some 
results of the tests of the various data windows are tab- 
ulated in Appendix B. 
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A block dlugrum for t ho bandpass SNR ostimation 
routine. SNRBPS, is t^ivon in Fipuros 10 and 11. Figure 
10 gives an overall view of the routine. Unless the group 
delay is known and is input to the routine, the first step 
required is to find an initial estimate for t^, denoted 
by T^. Then x(t) is delayed and blocks of data from 
X (t) and y(t) are transformed, windowing if specified. 

The power in x, is computed as shown, averaging over 

all the blocks of data processed. The total power, P^, 
is twice the value of P, defined by (48), and the power 
in y. P . is similarly found. Using delay values in the 
vicinity of computations of ( 2 R 2 > normalized by 

/FT- are performed for the number of Iterations specified 
by input parameter ILROAM. The maximum value obtained 
in these computations is found, and RTMAX is defined as 
twice R 2 fo simplify the notation. Finally, computation^ 
are done for estimating the SNR and for ERRN , the norma- 
lized mean-square error. All the final computations 
except the FFT are done using double precision arithmetic. 

Figure 11 provides a view of the method used to 
form the initial group delay estimate, The process 

involves first finding the peak value of the cross- 
correlation function R^y(T',. This process is identical 
to that described in the discussion of the lowpass SNR 
estimator. Once the delay value associated with this 
peak. is found, a block of data centered about 

is interval the Hilbert 
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t lu'n usi'd as a starting point f*'r i lu' dolay iti-rations 
porformod in t ho final <*ompnt at t I'tis . A listing t't t !i»' 
oomputor suluMUttno SNRIU'S is givon in Appondix 0 . 

All the considerations mentioned in the discussion 
of the Icwpass SNR estimation routine concerning the 
total number of data points used and the chviloe »'f the 
FFT block length apply as well to the routine for SNR 
estimation in simulations of bandpass systems. Larger 
sots of data generally result in mv're acourati' I'stimates 
for the SNR. Errors within llu' FFT increasi' as bU>ck 
si^.e increase.s. However, the use t'f v«'ry short bKu’ks 
can cause prt'bli'ius with t lu' slat iv'iiarity assunipt it'iis 
w*ith tilt' initial group delay t'stimatt'. Rt'sults t' f tests 
using varit>us blt'i'k It'iigtlis art' givt'ii in t lu' st'iM ivui I'li 
app 1 i ca t t t>ns 1 1 ' anaU*g systt'ins. 
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VI. .APPLICATION TO ANALOG SYSTEMS 

In order to study the aoouracy, rt>l i ab i 1 i t y , and 
range of the SNR estimation routines, eomputi'r simulations 
were performed for a number of analog .syst<'ms. Initially, 
test signals were combined with filtered noise, and the 
SNR estimation routines were applied to both lowpass and 
bandpass systems. Additive bandlimited white Gaussian 
noise was employed, and delays were prograimned into the 
data arrays in many cases to check the ability of the 
routines to correctly estimate the values of these 
delays. The test signals used included simple sinusoidal 
signals and iiK>re complex waveforms with Gaussian distrib- 
uted amplitudes. Tests were then made on a more compli- 
cated simulation of an FM communication system with 
phase-locked loop demodulation. 

In all the simulations performed for both analog 
and digital communication systems, filter models were re- 
quired. Therefore, a general lowpass and bandpass Butter- 
worth digital filter with variable order and critical 
frequencies was programmed for use in the simulations 
which follow. Some preliminary tests of the digital 
filter, as described in the following two paragraphs, 
were necessary before all of the results could be 
properly evaluated. 

To accurately determine the theoretical values of 
the SNR in the simulations, it is sometimes necessary to 
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compute the noise power at the output of a digital filter 
using measurements taken at the input. In order to per- 
form this computation, the noise-equivalent bandwidth of 
the particular filter involved must be found. In the 
simulations of digital communication systems discussed 
later, the error probability estimates can be properly 
formed from the SNR estimates only if the noise-equivalent 
bandwidths of the filters can be accurately determined. 
Therefore, tests were performed to measure the noise- 
equivalent bandwidth of the general filter for a range of 
filter parameters. 

For many analog filters the noise-equivalent band- 
widths are well-known, and tables of these values may be 
used to estimate the noise-equivalent bandwidth of the 
digital filter. However, in many cases, especially where 
a critical frequency of the filter exceeds ten percent 
of the sampling frequency, the digital filter noise- 
equivalent bandwidth departs significantly from that of 
the corresponding analog filter. This is particularly 
true when the filter order is low. The noise-equivalent 
bandwidth, even for a given digital filter, will vary 
somewhat depending upon the details of the implementa- 
tion used. Different software packages may compute the 
filter coefficients in a variety of ways, resulting in 
different errors. In addition, the different register 


lengths of various machines will also create some v'aria- 
tion in the noise-equivalent bandwidth from computer to 
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computer. Thus, in order to determine very accurately 
the noise-equivalent bandwidth for a tjiven digital tilter 
implemented on a particular machine, the bandwidth should 
be e.xper iment al 1 y measurt'd on the comput r whert' it is 
to be used. The noise-equivalent bandwidtij of the general 
filter used in all the simulations perfi>rmed on the 
IBM 370 was tabulated for a wide range of filter param- 
eters using a technique for experimentally determining 
the bandwidth. Appendix D explains the method used to 
perform the necessary calculations. A brief table is 
included to give an indication of the results which were 
obtained. A complete tabulation of the rather extensive 
results is not given, since it would not necessarily be 
applicable to a different software implementation, or 
even to the same software used on another computer with 
a different register length. 

In order to test the lowpass and bandpass SNR estima- 
tors in estimations of a simple analog system, various delay 
values were programmed into test signals. Then, filtered 
Gaussian noist' was added to these delayed signals. Both 
lowpass and bandpass systems were simulated, and the two 
routines for SNR estimation were applied to measure the 
SNR and t.. ’ ly values. The signals themselves were not 

passed thri>ugh filters in these tests because unknown 
delay and distortion would be imposed by the filters. 

This might yield interesting results, but it would makt' 
evaluation of the SNR estimation routines inpc'ssible 
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because the theoretical values for the delay and the SNR 


would then be unknown. These initial tests were intended 


to check the accuracy and the useful operating range 


for the estimators. The tests were not intended to 


measure the delays and accuracy of the filter models. 


A block diagram of the configuration used to perform 


these simple simulations is shown in Figure 12. The test 


signals used included waveforms with Gaussian amplitude 


distributions having lowpass and bandpass spectra for use 


nth the corresponding SNR estimation routines. Sinusoidal 


signals were also used with the lowpass routine. Sinusoi- 


dal modulation of a sinusoid was included in the tests of 


the bandpass routine. The additive Gaussian noise was 


filtered so that the noise bandwidth was somewhat greater 


than that of the signal. Lowpass and bandpass noise 


spectra were usually chosen to correspond with the type 


of estimator being used, although this is certainly not 


a necessary requirement. From Figure 12 it is clear that 


the signal component of the measurement data is identically 


a delayed version of the signal source. Thus, even though 


the tesr, signal is bandpass, subroutine SNRMSE may be 


used because the group and phase delays are equal. The 


delayed signals were generated in some cases by time- 


shifting either the signal or, in the bandpass case, the 


modulation. In other cases, such as with sinusoids, 


delays were easily programmed into the signal generator 


itself 
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Figure 12. Configuration for Filtered Noise Tests 
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The results of these tests are hiven In Table 111 
and Table IV. I. should be noted that the SNR results 
are accurate to within one dll in all cases. For the 
SNBMSE tests the delay estimates are all correct. He- 
cause srouv delay estimates SNltWS are determined from 

the magnitude of the complex envelope, delay estimates 
may appear to correspond to either positive or nepaiive 
values of gain. However, a negative gain may alternately 
be considered as positive gain with a phase delay of 180 
degrees. For . signal such as the sinusoidally modulated 
test signal, theoretically correct delay values occur 
every half cycle of the modulation. There were only two 
errors in group delay estimation, and clearly the SNB 

estimates were not seriously degraded. 

In order to test the lowpass estimator further and 
to apply the routine to a simulation of a more complex 
analog system, a computer simulation was performed tor 
an FM communication system. The simulation mode 
quired included an FM modulator model and phase-locked 
loop (PLL) models tor demodulation, in addition to the 
filter package used in the previous tests. Two PLL models 
were utilised, one of which incorporated a hard-limiter 
at the input, and one of which did not. Transient re- 
sponse tests and phase error variance tests were performed 
to verity the accuracy of both PLL models. Thest result, 
which are not reproduced here, indicated satisfactory 
of the mode Is i 66 i . 


performance 


TABLE III 

IMULATION TEST RESULTS WITH SNRMSE 




101 


r^ 

CN CM CM 


o 


I 


p: 

< 

GO 


CO — 

XJ 

w 

z 

CO 


JC 

-p 

H bc 
Cm C 
Cm 0) 


t> 

CM ^ ^ CM 

o o o o 

OOOCOOO«-«r-.r-.f-> 




CiCiCiCiOCsOOOOO 

r-tCMCOinrHCOmCO^lO^ 


ocooooooooo 

CMCO’VCOCM^QCOTfinO 


Oi 

CO 

CO 


• 

» 

• 

“ 

0 

X* 

Oi 


CM 

CO 


CO 

0 

in 

0 

in 

CM 

CO 

in 

CO 


§ ^ ^ 15 ; ^ ^ 00 00 00 00 


X 

0) 

p a 


CO 


CO 


p 

X c 

0 ) tx 

CO 


^^00 00 '^X*X *00 00 00 00 

r-tr-^CMCM-- — 

CDCDCMCMCDCCCCCMCMCMCM 


CO 

00 

00 

QO 

a 

X* 



Q 

0 

0 

0 


CM 

CM 

CM 

p-H 

f-H 




0 

0 

C^‘ 

0 


0 

0 

CD 


c 

CD 

1 


^ -2 ^ 
< < < < 
G Q 


S 


. , O C Q o O 


Ch 

0) 


D 

a 

B 

O 

o 


CO CO CO CO ^ 

cx:cx:a:cca: 


s s 

SH 

2 : 

< c 

a: o: 


< 

a 


< 


m; 

CO 

;p 


CO 

;-:) 


o 

CO 

c:2 


O 

CO 


X 

t 


1 c 


D4 C 

c 


V.*' 

>> 


a 

X 

c -- 

D w' 


”0 



XJ 

C P 

c 





Mi 

vTj S-. 

U 

^ c 


w w 

P 

X 


p r4 


CJ 

r-t 

sD 

r— 4 

u c 

< 

M. 0 


0 0 


0 CM 

X 

X 

X P 
-p 0 

G> 


u 

0 C 

0 

0 r-* 

u 

0 a 

c. 

CM P 


p 

c 

p p 


0 :3 

■f* 

E 

X 

0) 

X 

rH >> 

D 

a c 

oi 

P fA 

C 

P w 


p 

X 

:3 -f 

X 

£ 

ci 

0 

iD« 

CD 

03 

c c 

C 







>. >! 


p p 


.f-t 

cz 

c c 

-> 

.r^ r-^ 


03 *03 

p 

c c 



c;: 

• M. ^ 

■»-4 

C w D 
ctf Cm 

P 



CO CO CO CO 


CO 


CO 


CO 


CO 


7 : ' 

c c/3 

3 X 3 T 3 

*»-t p 

X O 


0000000000 o u o 

ssssssssssg^a^s^g, 

ilill§iiiiSi-B 5 i""" 

^ ^ HH ;:d P P P P 


r - 4 ;j) *— » 


P X 

X X 


X X 
:3 z 5 


0) 

p 


c 
ci *'-• 
a c/: 

o c 


uO 


:3 

K 


CMCO^iOCOC^OOO^O 


« 

» 

I CM 


* 

» 




C 

a: 


a; C; C3 c 

M. D 3 3 

rt a a 

* ♦ 

^ ♦ 

Ji- 









..,.,1 ■ 




.4 


102 


o 

t^ 

CO 

§ 


CQ 

S; 

A 

« 

2: 

CO 


> S33 

H-i H 


iJi 

03 

< 

H 


-3 

:=i 

CO 

M 

« 

CO 

pq 

Z 

O 

l-H 

H 

:=» 


CO 


e 

H 

CO CO CO CO CO TT CO 

COCOCOCOCDCCCOCD 

CD CQ r~^ r—t r-H r-H Tj^ CD O CQ 


COCOCOCOCOCOCOCO 


COCDCQCQCQCQCDCO 

H 

O iP 1 — ♦rHf-HiPr-*r— ( iPOOOOOO 

ft; 

rHO3T|*iO00C0a3^C3CMr-ICjTrTj*00 

►** 

< 2: 

OOOCCSiCiOilPOfPOli— iOiPQ3 

CO 

CiOfHCO Ci r-»ip caco^^ 



OQ 


•a 




Art 

C^CJ^^03CJJCOtnOl0005i-lCM^ 

2; 

00000>0>0000000i-<0j 

CQ 

W CJ r-( CO «(Hr-lr-( N CO ^ Tf 

x: 


-p 


H tkfi 

CQCD00 0003C^CDCQCQCDCOCOCOCOCD 


O>Orr^C3^C^0iOio>CiO>(j30)abO 

P^ a; 

0 

0 

0 

0 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 



CO 


r-H ^ 

oooooooo oooooo 

Oi rH 

ooooooooooooooo 

a 

ooooooooooooooo 

0 e 


H c^ 

ootnmoomoooooooo 

CO 

rH r-t T— < Cq rH OJ rH P nH rH rH rH 

p 


o 


«H iH 


w'^ 

hD i-q 1-3 -q 

c c 

<: <; c c <3 

O bO 

QQQQQC. jSSSSSSSS 

■H *H 

I— ^OOOOOOOO 

P CO 

OOOOOOQpPQppQpQ 

cd 

cococncQCQcn^zz^zzzSz; 

iH P 

opjDtDC3t3<;<<<<:<<<< 

^ CO 

^zzp2;;;z;o::a:cx;(i;tfp(:q:B:oc: 

T3 <U 

IP IP PI h-l ►H H-t 

0 H 

CO CO CO CO CO CO 

2 


C 


:3 

* 

(X 

rHCaCOTj<iOCOt^OOC3>OPCS|CO-^iO 


p p p p p p 


o 

4J 

<u 

Z5 

•o 


- 

cS 

cu 


u 

CO 

0 T3 


p 

(1) 

(L 

u 

p 

(I) 

a; 

sz 

p 

H 

^ rt 

o E 

O p 

• 

p p 

CQ 

CO 

a> CM 0) 

P o 

cu 

Pi 

I 

02; 
P CO 

OQ 

a 

p -c 

o 

p p 

o 

p p 

p 

o <a 


£ > 

CQ 


p 

+ |d 

c 

CO 0 

0 

CO 0 

p 

P c 

p 

W 0) 


x; 

p 

p 

:3 

O T3 

TD 

C 

0 

o; as 

s 

P 

cd CS 

p 

£ p: 

ctf 

p 


p p 

p 

CO O 

o 

0) 

CQ 

0) 

p 

>i ;3 

c 

cd P 

p 

P 

vO 

<D > 
TD 

(D 

P 

x: 

u o 

p 

P Q) 


O ^ 

p 


0 

bC O 


o 

Q) 

cd 

rP 

>> 

o 

p 

>» 

>>p 

o 

p cd 
P o 

p 

O P 

p 

p p 

ctf 0) 

x: 

o ^ 
P o 

(U 

P <D 

c 

<D x: 

o 

a p 

* 



5 




















103 


In the FM communication system simulation, several 
different combinations of system parameters were used. 
Simulations were performed using 10 Hz sinusoidal modula- 
tion and a 200 Hz carrier with deviation ratios of one 
and five. Actually, the frequency scale utilized was 
arbitrary. For example, if 10 MHz modulation and a 
200 MHz carrier were used, the simulation would remain 
the same except for scale factors. The frequencies used 
were chosen mainly for convenience in determining various 
system parameters. Two different sets of seeds were used 
for the noise generator. Then, a sample function from 
a lowpass Gaussian random process, ^ 3^13 ® 30 Hz, was used 
for the modulating signal with a deviation ratio of one. 
Once again two different seed sets were utilized for 
generating noise. Every simulation was performed using 
both of the PLL models. In each case the system was 
simulated for a total of one second with a sampling fre- 
quency of 10 KHz. Data were collected for the last 0.75 
second. Therefore, results are based upon 7500 samples 
in each case. A block diagram of the system under test 
is shown with the SNR estimator SNRMSE in Figure 13. 

The simulations were performed on two computer 
systems, an IBM 370 system at the University of Missouri 
and a UNIVAC 1108/1110 system at NASA-Johnson Space 
Center. All of the cases described above were tested 
using the IBM 370 system, and some of the tests were 
duplicated on the UNIVAC 1108/1110 system. Specifically, 
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tho simulat ti>ns usin>: sinusoidal imululation and a lit'via- 
t ii'n ratio of fjvo wort' pt'rfornu'd on both mai-iuiu's 
iH't h ril nioiit'ls. Tt'sts wt'i’t' niadt' usin>; i’ll natural fr»'- 
tiut'noit's ol f^^ - ut ;uui f^^ - 10(>. 

In all tt'sts. th»' lowpass SNH t i nia t i on j’outino. 
bN'RMSi:. was ust'd to nu'asurt' t ht' s i pna 1 - t o -no i s,' ratios 
at tho input to t ho Pl,l, dt'int>du lat or . vbNH)^. and at t lu' 
output of tho p»>stdoloot ion flltor. Altln>u>:h 

the signa.1 at the input to tho Pl,l< was bandpass, it 
may be seen from Figure 13 that there is no delay involved 
between the reference and the measurement data. Therefore, 
the routine developed for lowpass systems may be applied. 
Predetection slgnal-to-noise ratio. (SNR),^. wa.s varied 
over a range of values to include the region (’f thro.sht>id. 
These values and the resulting values of postdot «'ot imi 
signal -to- no iso ratio. (SNR)j^. woro tabulatod. This 
data appears in Appendix E. 

Tho thoorotioal value of whioh results from 

a given value of iSNRl.^. for an FM system using a oonven- 
tional discriminator with sinusoiilal modulatit'n is given 
by ib71* 


...1 


(SNR) 


P 


2W 


I 


2WP, 


(SNR). 


;ir 


er fe| /( SNR ' *■ 


8\V 


(SNin. 


This expression differs from that given by Ziemer and 
Tranter because a noise term is incliuU'd iu>r«' 
negl.'cted in t lu' refer,>nee to v'btain a simple 
mat ion. There are also notalional changes 


which is 
approx 1 - 
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where 

W = Postdetection filtt>r noise-equivalent bandwidth 
= Predeteotion filter noise-equivalent bandwidth 

and 

f^ = Peak frequency d('viation. 

For the case of Gaussian modulation, the corresponding; 
expression is [68]* 


W‘ 


(SNR)jj»| 2 2 ’ (SNRU 




u IB 


erfc[/(SNR)^] 


2/2W ^ 


-(SNR), 


-1 


2 

where o represents the variance of the modulating signal, 
m 

Theoretically, the PLL and the conventional discrimina- 
tor performances should be essentially the same for high 
values of the SNR. As the SNR is reduced, however, the 
PLL should extend to a lower value of (SNR),^ before 
exhibiting a threshold. Exactly this behavior can be 
seen in the curves of (SNR)^ versus (SNR),p which were 
plotted in Figure 14 through Figure 21. The threshold 
extension is more pronounced at the higher deviation ratio, 
just as it should be. This may be seen by comparing 
Figures 14 through 17 (deviation ratio = 5) with Figures 
18 and 19 (deviation ratio = 1). Figures 14 through 17 
also show the effect of decreasing the loop natural 

*This expression includes a noise term which is neglected to 
obtain a simple approximation '^y Taub and Schilling. Again 
there are notational changes. \ 













Figure 15. (SNR)^ Versus (SNR)^^, Case 2 
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frequency, As fj^ is reduced, the passband of the 

PLL is reduced. Less noise power is passed by the loop, 

and the threshold is extended more. This may be seen by 

comparing Figures 14 and 15 (f^^ = 100) with Figures 16 

and 17 (f = 40). Continuing to reduce the value of f . 

n " 

however, would eventually prevent the PLL from tracking 
the signal accurately because of insufficient demodulator 
bandwidth . 

The results obtained using Gaussian modulation are 
shown in Figures 20 and 21. The measurements of (SNR)^ 
appear to be low for large values of (SNR),p. This reduc- 
tion in accuracy was partially caused by distortion of 
the lowpass random signal in the postdetection filter, 
since a small part of the signal power fell outside the 
passband of this filter. There were also some problems 
with obtaining stationary signals due to the very low 
frequency components in the modulation. Even so, the 
results appear to be quite good below (SNR),p = 20 dB. 

The accuracy of the simulation should be improved by 
increasing the order of the filter used to generate the 
lowpass modulating signal, so that its spectrum would ro"l 
off more sharply. Reducing the amplitude of the lowest 
frequency components of the modulation should reduce the 
problems with stationarity and improve the results, also. 

The simulation results are in good agreement with 
the theoretical results for a true FM system. It is 
important to note that subroutine SNRMSE has provided 


accurate SNR estimates for valu€>s of the SNR from below 
0 dB to over 4 5 dB. Tlu'st' results indicat*' that the SNR 
estimation routine is valid and accurate. Tht'y alst* 
demonstrat*' the successful application of the SNR estima 
tor to a simulation of an analog? communication system. 
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VII. APPLICATION TO DIGITAL SYSTEMS 

One of the major reast>ns for (ievelopin^^ the SNR 
estimators was to facilitate the evaluation of di^Mtal 
communication systems usinp computer simulation. Thus, 
to test the routines, a number of digital systems were 
simulated, and their performances were evaluated using 
SNR estimation. These results were then compared with the 
theoretical performances. The systems simulated include 
coherent ASK and FSK systems and a noncoherent FSK system, 
all with simple additive Gaussian noise channels. Addi 
tionally, a noncoherent FSK system was tested in a Ray- 
leigh fading environment, and in another simulation 
diversity transmission was utilized to improve performance 
in the fading channel. In all cases the systems being 
simulated were chosen so that the mapping from the SNR 
to Pj. could be analytically determined. This was done 
so that could be estimated using the SNR estimate for 
comparison with a direct count of the errors produced 
by the system. 

In order to provide a somewhat realistic simulation 
of the actual system in each case, additive white Gaussian 
noise was combined with the signal in the channel, and 
the signal-plus-noise was filtered by a predetection 
filter. This procedure allowed unknown group and phase 
delays to be introduced by the predetection filter, and 

routine for bandpass systems was 


the SNR estimation 
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required. Some early attempts to use the lowpass ver- 
sion of the estimator showed SNRMSK t be completely in- 
adequate with such systems as, of iu>urse, would bt* 
expected. 

There is one theoretical problem that results wht'n 
filtering the signals in this manner which should be 
mentioned. The theoretical mapping from the SNR to 
is generally developed under the assumption that the 
only errors present in the noisy waveform are those pro- 
duced by the additive random noise. The signal is usually 
assumed to pass through the predetection filter without 
distortion. In practice, the filter generally distorts 
the signal to some degree. The amplitude response of the 
filter is not perfectly flat over the passband, and the 
phase characteristic is not perfectly linear. Additionally, 
there are usually some spectral components of any prac- 
tical signal which fall outside the passband, and inter- 
symbol interference can become significant. In order to 
justify neglecting errors in the mapping from the SNR to 

P„ due to these effects, a simulation was performeu for 
h 

each system with the Gaussian noise set to zero. The 
SNR estimator measured errors in the system due to all 
sources other than the Gaussian noise. It was determined 
in each case that the Gaussian noise error would exceed 
all other errors by at least ten dB (and usually much more), 
even for the lowest values of Gaussian noise power to be 
considered. Hence, the assumption that the total noise 
power is essentially the result of only the additive 
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Gaussian noise is justified, and the theoretical mappinfj 
from the SNR to P„ is applicable. To further substantiate 
these arguments, additional simulations were performed 
for the coherent ASK and FSK systt'ms in which the Gaussian 
noise was first band 1 ini i tt‘d by a filter identieal to tin* 
predetection filter befort* adding it to tlit* signal. Th<Mi 
the predetection filter was removed from the signal path 
altogether. The results of these simulations are tabulated 
and plotted with the other results for the coherent ASK 
and FSK systems. It can be seen that there is no notice- 
able difference in the accuracy obtained with the two con- 
figurations . 

Block diagrams of the simulations used for the coherent 
ASK system are shown in Figures 22 and 23. The modulation 
source utilized to provide input symbols to the binary 
ASK modulator was a PN sequence generator producing pseudo- 
noise symbols in a sequence with a period of 1023 symbols. 
The ASK modulator performed on-off keying of a carrier. 

In one case, white Gaussian noise was added, and the 
noisy waveform was passed through a bandpass predetection 
filter. In the other case, the same filter was applied 
to the noise before it was added to the signal , and the 
signal was not passed through a filter. In both cases, 
the output was provided to a coherent ASK demodulator and 
to a PLL which tracked the carrier component and provided 
a phase reference for the coherent demodulator. The de- 
modulated symbols were compared with the original modula- 
tion symbols, and the errors which occurred in operating 









Figure 23. Simulation of the Coherent ASK System Without a Prede 
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th.- system were counted. Dividing the number of errors 
by the total number of symbols proi-.-ss.-d yi.'lds an 
estimate for the probability of .ri-i.r obtained by a 


airoot count , Pj^^.C count ) . 


Data were collected at the input to the di'niodii 1 at.>r 
where an estimate of the SNK was desired. Data were also 
collected at the modulator output to serve as the reference 
signal. Subroutine SNRBPS formed the SNR estimate for the 
waveform at the demodulator input. From this estimate, SNR. 
an estimate can be obtained for the ratio of the average 


energy per transmitted symbol. E^. to the noise power spec- 


tral density. N^. Then the well known mapping from E^/N^ 


to P for a coherent binary ASK system operating in addl 
£ 

tive Gaussian noise, [691 


Tj, = 5 erf 


to A 


s 


2 N 


(60) 


o 


where 


erf c(x) 


2 _ 


e cia 


can be applied. 


In order to find E /N^ from the SNR, denote the 


time duration of an information symbol by and the 
average power present in the signal at the point of mea- 


surement by P . Then the total average energy per symbol 


is given by 




t 

I 


, ..(Hr M ' ■-* ■ i- 
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E = P T . 
s s s 


Denote the noise-equivalent bandwidth ot the filter bv 

B Then the total noist' power prt'st'nt at the point 
n 

of measurement is given by • Thus. tht> SNR at the 

point of measurement can be written 


SNR = 


o n 


£ 

s 


o s n 


Therefore 



^ ■ <SNR)(T^B^> 
o 

Hence, the ratio E^/N is related to the SNR by a pro- 
portionality constant which is the time-bandwidth product 
of the system. The relation for P^. becomes 

P^ = I erfe /|(SNR)(TgB^) . ( 61 ) 

Inherent in the theoretical development of this map- 
ping is the assumption of a perfect phase reference in 
the coherent detector. Although the phase reference pro- 
vided by the PEL is not perfect, the phase error variance 
was measured for each simulation and was found to remain 
below 0.03 radian^ throughout. Hence phase errors may 
reasonably be assumed to be quite small, and the mapping 
based upon a perfect phase reference should yield a very 
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dost' ayiprox i mat i on to tho oorrt'ct mapping: tor t lu* systt'm 
being simulated. 

Simulations wt're pertormed tor this system iu<>r a 
rangt' t SNR values wlu're a direet error eount was tt'asible 
The various system parameters utilized an' tabulated with 
th«' results ol the .simulatii>ns in A{>pomlix K. A plot ot 
the results is shown in Figure 21. Points are plottt'd 
showing the relationship between the estimated SNU and 
the error probability obtained by a direct error count. 
Every data point represents the results of a different 
simulation. For each point, the value of the abcissa 
is the SNR estimate scaled by the time-bandwidth product, 

A 

(SNR)(TgB^), expressed in dB, and the value of the ordinate 
is the probability of error computed by a direct error 
count for that same simulation. The number of symbt>l 
errors that occurred in each case is shown next to the 
point in order to reflect the reliability of the P 

E 

estimate obtained by direct error counting. The curve 
showing the theoretical mapping from the SNR to Pj. is 
plotted for comparison. 

The distance along the ordinate between each point 

and the curve is the error in estimating the probability 

of error for that simulation that would result if P were 

E 

determined by estimating the SNR and mapping this to an 
estimate for Pg. Inspection of the plot shows good agree- 
ment between the Pg estimates which are found using the 



Probability of Symbol Error 
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10.LOG^^j[ (SNRXT^B^) 1 



Figure 24. Pj^ Versus SNR, Coherent ASK 
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SNR fstimation routines and the results obtained by dirt».t 
error counting. In all cases where more than five errors 
were made, so that the error count may Ih> considt'red im'an- 
ingful. the values obtained for the probability of error 
by counting and through SNR estimation differ by less 
than a factor of l.S. This level of accuracy is more 
than adequate in most applications where an estimate ol 


pjg is required. 

An important point about the data required to produce 
the estimate for Pg using the two techniques should be 
made here. Each ASK simulation processed a total of 30.700 
samples to simulate the 1023 information symbols used. All 
these data were used in finding the probability of erro»* 
by direct error count. The SNR estimate was produced 
using only 20,000 of these samples. Applying the mapping 
from the SNR to Pg for the system provided an estimate for 
Pg. Thus the simulation time required for this case cunild 
be reduced by more than one third with no change whatsoever 
in the P, estimates obtained through SNR estimation. Fur- 
ther, all the tests which have been performed to check 
the accuracy of the SNR estimator indicate that SNR estim- 
ates could have been obtained with only 10.000 samples, 
or even less, which would probably differ from the 
estimates actually obtained by no more than a few tenths 
of a dB. It seems likely that good Pg estimates for this 
system could have been obtained through SNR estimation 
for one third or less the cost in computer time which was 
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used. Similar results should be possible tor other sys- 
tems which can be accurately simulated, providing the 
mapping from the SNR to Pp is known. In the cases of 
the lower error probabilities obtained, such a reduction 
in simulation time would provide an insufficient numbt'r of 
errors to yield a reliable e.stimate of the error prob- 
ability using a direct terror count. This sort of reason- 
ing was experimentally verified in the FSK simulations 
which are discussed next. The estimates obtained 
through SNR estimation in those (FSK) simulations are 
formed using well under one half the data used for the 
direct error counts. The agreement between the results 
and the theory is even better for the FSK simulations than 
it is for the ASK case. 

The importance of these arguments is obvious when 

the extension to estimates of P„ in ranges several orders 

t, 

of magnitude smaller are considered. For the same number 
of samples, the reliability of the SNR estimate is 
almost unchanged as the SNR is increased approximately' 
ten dB in most systems to obtain this extension. This 
is seen in the tests of the estimators and in the confi- 
dence intervals which were derived in earlier sections. 

For these higher values of the SNR, a given error in 
•SNR generates an error in the estimate for P„ larger than 
at low values of the SNR. A brief study of this effect 
is given in Appendix F. The filtered noise test results 
obtained in the previous section show that the SNR 
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estimates are usually accurate to within a few tenths of 
one dB. In most digital systems, these measurements would 

lead to Pt^ estimates well within one order of magnitude 

E 

of the theoretical value. This accuracy is sufficiently 
adequate to provide a useful estimate of in most appli- 
cations. The computer time required to estimate P^. using 
direct error counts in this range usually exceeds prac- 
tical values by a few orders of magnitude, since, for 
a given value of the expected number of errors 
generated is directly proportional to the simulation time 
used. 

Thus, the SNR estimation technique for determining 
Pg is directly verified, at least over the range of 
shown in the plots. There appears to be no reason why 
the method should not also give good results for the 
lower ranges of Pj, typical in most practical systems. 

A similar simulation scheme was used to model a 
binary coherent FSK communication system. Block diagrams 
for the FSK system are given in Figures 25 and 26. The 
configurations are similar to those for the ASK system 
except for minor changes. The ASK modulator and demodula- 
tor are replaced by an FSK modulator and demodulator . 

Also, two phase-locked loops are required instead of 

one to supply phase references to the coherent demodulator. 

Again, a PN sequence is used as a modulation source. 

A single predetection filter is used with sufficient 
bandwidth to cover both frequency channels. The use of 















Figure 26. Simulation of the Coherent FSK System Without a Predetection Filter 











two filters with narrower bandwidths, one for each fre- 
quency channel, would yield a system with performance 
superior to that of the one simulated, since the SNR for 
each channel would be improved. However, this change 
would have no significant effect on the validity of the 
tests of the SNR estimator. It would introduce the prob- 
lem of matching the possibly different delays in the two 
filters to obtain the synchronization necessary at the 
end of each symbol time when a decision must be made 
about which symbol was sent. Also, the equivalent con- 
figuration which avoids filtering the signal, as shown 
in Figure 26, could not be so directly simulated using 
the version with two filters. Therefore, for convenience, 
the single filter was chosen. 

Simulations were performed over a range of SNR values 
as was done for the ASK case. System parameters and re- 
sults are given in Appendix E. A plot is drawn just as 
before to show the .elationship of the SNR estimates, 
the probability of error determined by direct error count, 
and the theoretical mapping between the SNR and P . The 
theoretical mapping is identical to that given by (61) 
for the coherent ASK system [10], These results, shown 
in Figure 27, indicate excellent accuracy in the SNR 
estimator. In this case the phase error variance re- 
mained below 0.04 radian^, indicating nearly perfect 
coherence. The values obtained by the two methods for P^, 
again differ by less than a factor of 1.3, except in cases 
where less than five errors were counted 
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Discussion of the accuracy of the methods developed 
to evaluate the performanct' of digital communication sys- 
tems is most directly acctnnp 1 ished by considering the 
error in estimating the probability of error. For 
example, in the ASK and FSK systems discussed, the varia- 
tion between the probability of error obtained by direct 
error counting and the Pg estimate obtained through SNR 
estimation was noted to be less than a factor of 1.3 
in all cases where more than five errors were counted. 

But the accuracy of the SNR estimator is roost easily 
described in terms of the error in the SNR estimate, 
not in terms of the error generated in the Pg estimate. 

It is useful to have a set of curves available which 
maps the error in the SNR estimate to the error in the 

P„ estimate. Such a mapping is obviously system depen- 
E 

dent and cannot be obtained in general. However, the 
case of coherent systems operating in additive Gaussian 
noise is quite common, and it is helpful to derive curves 
for these systems. In Appendix F curv'es are obtained 
which cover some of these systems, and a convenient way 
of expressing the error in estimating Pg is suggested. 

The data obtained from the coherent ASK and FSK simula- 
tions is plotted with the curves to show the area in 
which these simulations were operating. 

The next system which was simulated is a noncoherent 
FSK system. The simulation is the same as that for the 
coherent FSK system except for the form of the demodulator 




noncoherent demodulator 



f the Noncoherent FSK System 
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Assuming the predeteotion filter passbiind includes th*‘ 
passbands of the individual filters which follow, as 
it must for pn>per operation, the p»> rf (uanance (>l l ht> 
system is essentially t lu> sanu' whetluT or not the j)re- 
detect ion filter is ustul. Without it. howevt'r. the signal- 
to-noise ratios would have ti> bt> measuri'd at the outputs 
of the demodulator bandpass filtt^rs. This prt's»»nts a 
problem in determining a r».' ft'renct* signal bt.*oaust' t'ach fil- 
ter passes frequencies corresponding to only one of the two 
information symbols. In order to measure the SNR at the 
output of one of the demodulator filters, a special ref- 
erence signal would have to be generated corresponding 
to the theoretical output of that filter. This might be 
possible, but the predetection filter approach appeared to 
be simpler. Since the quantity to be determined is 
E /N , the calculation is performed just as in every other 
case, multiplying the SNR by the time-bandwidth of the 
filter involvf^d, in this ease tht' prt'det t ion filter. 

Once is found, the expression [71] 


E 

s 



can be used to estimate 

E 

The results of simulations performed for the system 

are tabulated in Appendix E, and results are illustrated 

in Figure 29 • Once again, the results appear to be in 

rather good agreement with theory, and the SNR estimator 

is shown to be capable of producing good estimate's fcir P« 

L 


Probability of Symbol Error 
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Values for the probability of error derived usinp; the 
two methods differ by less than a faetor of 1.5 in all 
cases for simulations of this system. 

In the ne.\t set of simulations, the sanu' noncoht'rent 
FSK system was used. liowt'ver. the channel model included 
a multiplier fiu- the signal which simulated Raylc'igh fad- 
ing. A block diagram for the system is shown in P’igure 
30. The tap for the reference signal tor SNRBPS is taken 
at the output of the random gain multiplier. The SNR 
value which is required to compute Pj. for the system is 
based upon the average symbol energy after the Rayleigh 
fading. If the reference for SNRBPS were taken at the 
output of the FSK modulator, all error produced by the 
variations in signal amplitude due to fading would be in- 
cluded as noise in the SNR estimate. 

The Rayleigh random variable used to simulate fading 
was generated as the square root of the sum of the squares 
of two independent, equal variance Gau.ssian random \ari- 
ables. In order to achieve an intuitively appealing model 
for the fading, it was desired that the fade durations be 
long enough so that the signal amplitude would remain more 
or less constant over the width of a single information 
symbol. Yet the duration was required to be sufficiently 
short that several fades would be included in a single 
FFT block size suitable for use in the SNR estimation 
routine. Otherwise, the st a t ionari ty assumption.s required 
of the signal for proper operation of SNRBPS would be 
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violated. The adjustment of fade duration was performed 
by controlling the spectral content of the indep(*nd(*nt 
Gaussian random variables utilized. Independent Gaussian, 
equal variance random variables were used as inputs to 
two identical bandpass filters. Since the two filters 
are linear systems, the outputs remain Gaussian. Thus 
two independent Gaussian random variables of equal variance 
were available at the filter outputs at each sample time 
to form a Rayleigh distributed random variable. The 
approxiotate duration of the fades was controlled by iden- 
tical adjustment of the two filters. The desired fade 
duration was achieved by experimenting with the filters 
and observing the results. Statistical tests were per- 
formed to verify the independence of the Gaussian random 
variables between the two filter outputs and to verify 
that the distribution of the final random variable 
generated was a good approximation to a Rayleigh distri- 
but ion . 

The value of P„ for the Rayleigh fading channel is 

tj 

given by [72] 



where E is the time-averaged symbol energy present at 
s 

the demodulator input and is the powe' spectral density 
of the white Gaussian noise. This mapp' ng was used to 
obtain estimates of Pr, from the SNR estimates. A large 
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number of fades were included in the data stream used to 
compute the SNR, so that an accurate value could be ob- 
tained for the time-averaged symbol energy. A table of 
the results of these simulations is given in Appt'ndix I). 

A plot similar to those described for the other sysftmis 
is given in Figure 31. The P^, values obtained by SNR 
estimation and by direct error count differ in all cases 
by less than a factor of 1.2. Once again the SNR esti- 
mation routine appears to be extremely accurate. 

In the final sirailations performed for digital 
coonunication systeiis, a diversity tran^lsslon system 
was simulated to combat the Rayleigh fading channel. It 
was assumed that a nunA>er of independent Rayleigh channels 
having Identical statistics were to be utilized. Spatial 
diversity was actually simulated, although the results 
would theoretically be the same if the system employed 
time or frequency diversity, as long as the assumptions 
of independence and identical statistics remain valid. 

The modulation scheme employed was noncoherent FSK. 

The system configuration for each spatial channel was 
identical to that of the single channel system just de- 
scribed in the previous simulations except for the decision 
mechanism in the demodulator. In order to decode a symbol 
in the diversity system demodulator, first the outputs 
of the two envelope detectors in each binary channel are 
squared. Then the squared outputs corresponding to the 
same symbol in each channel are summed over all the channels. 


10-LOG^Q[(SNR)(TgB^) J 

t 6 8 10 12 14 16 18 

f ] I 1 " » I ' ■ 1 1 ^ ’ * * 

NONCOHERENT FSK SYSTEM 
WTTH RAYLEIGH FADING 

SIMULATION RESULTS 

theoretical Pj, versus SNR 

NOTE: (SNR)(T B^) = 


o224 


Figure 31. Pp Versus SNR, Rayleigh Fading Channel 
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Finally the larger of these two sums is chosen to define 
the decoded output symbol. The structure of t h i .s system 
is illustrated in Figure 32. 

The effects of diversity transmission can be .seen 
by assuming a given amt>unt of t'nergy is availabU' for 
eac ii modulation symbol . Holding this total .symbol en<'rgy 
constant, system performance is ob.served a.s tlu? t'nt'rgy 
is divided among an increasing number of independent 
channels. This procedure was followed in the simulations 
performed. The value of Pg in this system is given by 
C73] 

Pe * (1-p)’^ 

j=0 J 


where 

1 

■ [''"oj 

L = the order of diversity 

and Eg and are as defined earlier, but remembering 
that Eg represents the total energy summed over all the 
channels . 

In these simulations the ratio E /N was to be fixed 

s o 

at 15 dB. The SNR was est imated for only one channel. 

It was assumed that the other channels were operating 
at the same SNR. Each order of diversity from one through 
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four was simulatod. In oarh oaso dat a woro t“(^l loctod from 

one channel to estimate the SNR, and a computation was 

performed to obtain E /N . The vaUu' was lu‘ld c h)s«‘ to 

s' o 

15 dU for each order of diversity by usinp the SNR 
estimator to check performance. The number of errors 
which occurred in each case was used to obtain the proba- 
bility of error. The results are tabulated in AppendLx 
E. A plot was constructed showing the value obtained 
for the probability of error at each order of diversity. 

A curve was drawn indicating the theoretical value of 
p£ for each order of diversity, assualng that la 

dB. Of course, the curve is only defined for a practical 
system at Integer values of diversity . The plot is given 
In Figure 33. The results again are quite good. For 
diversities of orders less than four, the values of P^. 
miss the theoretical curve by less than a factor of 1.1. 
In the fourth order diversity system, the factor is about 
1.7. 

In this case the SNR estimation routine was not 
utilized in order to measure an unknown but fixed simula- 
tion. Rather, the routine was used to measure the SNR 
of a simulation to maintain the desired SNR specification 
for the system. Adjustments of the noise were checked 
using SNRBPS until the desired SNR was achieved. The 
SNR estimation routine appears to have been accurately 
monitoring the SNR, since the number of errors that 
occurred are in good agreement with the theory. Thi.s 
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particular application demonstrates another way the SNR 
estimation routines may be helpful in evaluating the 
performance of a communication system tlirough computer 
simulation. 


VIII . 


CONCLUSION 


As diijital computer simulation has be^un to play a 
major role in the analysis of many systems, it has been 
seen that a general criterion for performance evaluation, 
which is suitable for use in digital simulation, is needed. 
In many cases a useful criterion is the mean-square error 
between a measured waveform and some reference. In other 
cases some other figure of merit for a system can be de- 
rived using the mean-square error. In conjunction with 


measurements of power, the mean-square error can be used 
to define signal-to-noise ratios for arbitrary waveforms, 
where the error is considered to be noise. Other param- 
eters of general interest, such as system gain and system 
delay, can be defined by using measurements of the mean- 
square error or the signal-to-noise ratio. 

In this work, methods were developed for estimating 
all these quantities in simulations of lowpass and bandpass 
systems. The technique requires that a reference waveform 
be available which differs from the ideal signal component 
of the waveform at the point of measurement only by an 


unknown scale factor and. in the lowpass case, by an un- 
known time delay. In the bandpass case, both group and 
phase delays can be arbitrary. The method can be applied 
to digital simulations of any wide-sense stationary system 
where this reference is available. The measurements of 
power, gain, delay, and signal-to-noise ratio find many 
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applications in simulations of both analog and digital 

communicat ion systems . 

In order to obtain these estimates, a m(>thod was 
developed to determine the amplitude scaling; and the 
group and phase delays that must be imposed on the ref- 
erence to form a signal estimate which minimizes the mean- 
square error with respect to the measurement waveform. 

Two software implementations were developed. One is for 
the case where group and phase delay are equal, which means 
the signal estimate may be formed by simply asqjlltude 
scaling and time-shifting the reference. This case is 
usually associated with lowpass linear systems. The other 
implementation, which allows for arbitrary group and phase 
delays, reflects the processing generally encountered In 
bandpass linear systems. Both of the SNR estimators were 
proven to be ML estimators in systems with additive white 
Gaussian noise. Curves showing confidence intervals were 
developed for the lowpass SNR estimator in this case. 

These curves indicate that excellent SNR estimates can 
usually be formed using data sets of a practical size. 

Development of confidence intervals for the bandpass 
estimator proved to be too complicated for simple mathe- 
matical analysis. Further work could be done in deriving 
such curves. It would be helpful to have a more thorough 
statistical study of the estimators. It may be possible 
to determine whether the estimators are ML for channel 
models other than the additive white Gaussian channel, and 
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confidence intervals would be useful for these channels, | 

also. I 

I 

The SNR estimators were first tested in applications 
to analog systems. Results obtained from simulations of 
a filtered noise system and a more complex FM communi- 
cation system indicate that accurat'- estimates can b<* 
obtained for a wide range of signal-to-noise ratios. In 
some cases, the SNR estimator is useful in checking the 
validity of the simulation Itself. For example, in the 
FI system which was simulated, the SNR measurements helped 
to verify the proper operation of the PLL d^iodulator 
simulation model. 

A survey of the techniques currently available for 
evaluating the performance of digital coimnunication sys- 
tems through the use of computer simulation revealed that 
the methods generally in use are not suitable for evalua- 
tion of most practical systems. The parameter most often 
desired in evaluating these systems is the probability of 
symbol error, P^.. The estimation of P^, based upon count- 
ing the errors which occur, or extensions of error counting 
techniques, require too much computer time to be feasible 
for values of Pg less than about 10 which are typical 
in most systems. 

In many cases a mapping from Pj, to the SNR at some 
point in the system can be found, or at least approximated, 
using either mathematical analysis or empirical data. In 
such systems. P^ may be estimated through the measurement 
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of the SNR. In some cases the SNR itself is a useful 
parameter for evaluating performance. 

The technique developed for SNR estimation was 
applied to simulations of a variety of digital communica- 
tion systems. Estimates of Pj. obtained through SNR mea- 
surements were compared with results obtained by direct 
error counting. The tests directly verified that accurate 
estimates for Pg can be formed using the SNR estimators 
over the range of the SNR associated with error probabili- 
ties where error counting is feasible. Since measurements 
can be accurately made for a much higher range of the SNR, 
all indications are that accurate estimates for Pg can be 
found over ranges of Pg where the error counting methods 

cannot be applied successfully. 

In short, the results of all tests show that the SNR 

estimation methods developed accomplish their task quite 
well. Accurate measurements have been obtained for a 
wide range of signal-to-noise ratios when applied to 
simulations of both digital and analog communication sys- 
tems. The accurate estimation of Pg has been directly 
verified for error probabilities greater than 10 ^ through 
simulation. In almost all cases the probability of error 
was estimated to within a factor of 1.5. It seems 
accuracy similar to this should be expected for the much 
smaller values of Pg common in most practical systems. 

The accuracy of the Pg estimate would be decreased slightly 
in many cases becuase. as the SNR increases, a given error 
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in the SNR estimate usually produces larger errors in the 
estimate for 

The method is economically ti^asible. Tht^ amount ot 

computer time required to dt'tt'rmiru' an SNR t^sl imate tor 

-6 

product nj^ a Pj^ estimate of 10 is t^ssentially thi^ same 
as for producing a Pj, estimate ot 10 , since tlu' riM^uiri'- 

ment on computer time is only a function of the numbt^r of 
samples being processed. The various measurements tor 
delay, gain, and correlation which are generated in the 
process are often helpful, also. 

The techniques developed seem to offer a versatile 
tCK>l for the accurate quantitative evaluation of many dif- 
ferent systems through the use of con^uter simulation. 



2 
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APPENDIX A 
SUBROUTINE SNRMSE 

Subroutine SNRMSE is a FORTRAN IV implementation of 
one of the SNR estimators developed for use in computer 
simulation. It is applicable to systems where the group 
and phase delays of the measurement data with respect to 
the reference data are equal. These are usually, but 
not necessarily, lowpass linear systems. For proper 
estimation of delay, the data should be wi de-sense sta- 
tionary. 

The measurement data set of K saoq;>les is input to 
the routine in the Y array, and the X array contains the 
reference data. Included la the argunent list are the 
various work arrays and dimension values required and 
input parameters to specify various processing options. 
These inputs include the choice of the FFT block size to 
be used, whether an initial delay estimate is to be deter- 
mined by the routine or specified in the argument list, 
the choice of the number of iterations to be performed in 
choosing the final delay estimate, whether to apply a data 
window, and whether the cross-correlation function is to be 
printed, plotted, both or neither. All the arguments of 
the subroutine are listed with explanations in the comments 
at the beginning of the program listing which is included. 

The execution of SNRMSE does not destroy any of the 
input data. Thus mvltiple calls of the routine or succes- 
sive processing of the data set by other routines is 
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possible. The routine provides a printout of estimates 
for power in X, power in Y, system gain, system delay, 
the peak of the cross-correlation function, normalized 
mean-square error, and the SNR estimate expressed in dB. 
Most of these values are also returned in the argument 
list so that they may be used in any successive process- 
ing. Also printed are a list of the values of the input 
parameters in effect upon execution of the routine, the 
normalized cross-correlation values used to select the 
final delay estimate, and the means and variances of the 
FFT data blocks ^icb were confuted in testing stationar- 
ity . 

If any stationarity tests were failed, a warning 
message is printed to that effect. The tests are made 
by computing the mean and variance of each FFT data block 
and comparing these values with those of the entire X 
or Y data set. If a value outside a certain experimentally 
chosen range is found, the warning is generated. To aid 
in locating the problem which generated this warning, 
immediately following the warning an output of the form 
IFLAG = I will be printed, where I is a signed integer. 

The magnitude of I indicates the number of the FFT mea- 
surement data block which caused the warning, except that 
a value of 10,000 denotes the reference block. A minus 
sign indicates a variance test was failed. Otherwise, 
a test of the mean was failed. If more than one test was 
failed, only the last failure v/hich was detected is 
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indicated by IFLAG. A warning statement does not neces- 
sarily mean that the results are in error, only that the 
stationarity is questionable. This often causes the 
delay estimate, and hence the other results, to be 
erroneous. If the correct delay value has been found in 
spite of the nonstat ionar i ty , then the results art' all 
valid. 

The subroutines CORR, GRAPH, and DATAPT, which are also 
listed, are called by SNRI4SE and require no attention by 
the user, except that they must be available to SNRUSE. 

The subroutine CORR is used to obtain the cross-correlation 
function, and GRAPH and DATAPT are used to graph this 

function if specified. Also, the FFT routine HARM must be 
available for use in SRRMSE. 

One parameter which needs some explanation is ISIGNR. 
This input specifies whether the cross-correlation peak 
to be found in determining delay is to be positive, nega- 
tive, or the largest absolute value computed. For most 
cases the user should specify ISIGNR = 0, indicating that 
the absolute value peak is to be used, and the routine will 
execute properly. Only for a special class of signals 
does the user need to specify whether a positive or nega- 
tive peak is required. If the signals being analyzed 
possess the proper symmetry such that the cross-correlation 
function theoretically has equal magnitude positive and 
negative peaks, then the user must specify which sign is 
to be used in order to obtain correct delay estimates. 
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This would be the case for signals with half-wave sym- 
metry. It should be noted that in such cases the esti- 
mates for the SNR, power, and error are valid regardless 
of the specification chosen for ISIGNR. Only if the user 
is concerned with delay or gain estimates in a system 
with theoretically equal magnitude positive and negative 
cross— correlat ion peaks, does ISIGNR need to be specified 
as 1 or 0. In such a case, a 1 should be used if the 
system gain is positive, and a -1 should be used if the 
system gain is negative. 
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APPENDIX B 

COMPARISON OF DATA WINDOWS FOR USE IN SNRBPS 

Initial testing of the Hamming data window for use 
in SNRBPS showed that a significant improvement in the 
accuracy of results could be realized compared to 
processing without applying the window function. 

These results prompted further testing of a variety of 
data windows. 

The tests were made using a double sideband signal , 
specifically a product of sinusoids. Amplitude scaling 
and delays were progranmed into the measurement signal, 
and additive Gaussian noise was applied. Tests were 
made for SNR values of 17,7 and 36.9 dB and for a very 
large value, which was too great to be measured accurately 
by the routine. Results are tabulated in Table V. Seven 
windows, in addition to rectangular, were tested. 

The results show all windows used work well for the 
lowest value of the SNR. At 36.9 dB, many window func- 
tions yield very poor results. Two Kaiser windows, the 
Hanning window, and the triangular window all remain 
within one dB of the correct SNR value. For the 
extremely high SNR test (on the order of 70 dB) , results 
are not tabulated, since no accurate measurements were 
obtained. However, the Hanning window provided an 
estimate in excess of 60 dB, which was more than 5 dB 
greater than any other window. On the basis of these 
results, the Hanning window was tentatively chosen for 


TABLE V 


RESULTS OF TESTS OF WINDOW FUNCTIONS 


Window SNR Estimates 


Function 

True SNR 
= 17.7 dB 

True SNR 
= 36.9 dB 

Rectangular 

16.4 

22.1 

Blackman 

17.4 

30.0 

Kaiser, 8 > ir 

17.5 

30.8 

Kaiser, 0 » 3w/2 

17.7 

36.3 

Kaiser, 8 * 2tt 

17.9 

37.6 

Haimning 

17.7 

34.8 

Hanning 

17.7 

37.6 

Triangular 

17.8 

37.4 
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use in SNRBPS. Subsequent testing verified the overall 
accuracy and reliability of estimates lormed using this 
window on many different signals. Therefore, the 
Hanning window was permanently installed in subroutine 
SNRBPS for the data window proce.s.si ng option. 
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APPENDIX C 
SUBROUTINE SNRBPS 

Subroutine SNRBPS is a FORTRAN IV implementation 
of the second SNR estimator developed for use in computer 
simulations of communication systems. It is applicable 
to systems where the group and phase delays of the mea- 
surement data with respect to the reference data are 
arbitrary. These arbitrary delays are associated with 
bau(4}ass systems. The routine will yield valid results 
only with signals which carry all information in the 
envelope of the waveform. In such cases only the envelope 
of the signal is required to be distortionless to retain 
all the information; delay of the carrier itself is un- 
important. The routine should be applied to systems 
having some carrier frequency displaced from zero. For 
proper estimation of group delay, the data should be wide- 
sense stationary. 

Input to the routine is essentially the same as for 
SNRMSE. The delay included in the argument list refers 
to group delay. There is an additional input parameter 
for SNRBPS which can be used to limit the maximum delay 
value to be considered in estimating group delay for the 
system. All the arguments for the subroutine are listed 
with explanations in the comments at the beginning of the 
program listing which is given. 

Like SNRMSE, the execution of SNRBPS does not destroy 
any input data, so successive processing of the data sets 


is possible. Outputs provided by SNRBPS are identical to 
those of SNRMSE with three exceptions. The delay value 
given denotes group delay. (No phase delay estimate is 
produced.) The peak of the envelope of the cro.ss-correl at ion 
function .s provided, instead of the peak of the 
cross-correlation function itself. Finally, magnitudes 
of normalized cross-spectral power estimates are printed 
instead of normalized cross-correlation values, since these 
power Beasuresents are used in estimating group delay. 

Stationarlty tests are applied and warnings generated 
exactly as in the lowpass routine. The subroutine CX)RR2 
is listed following SNRM»S. It is called by SNRBPS to 
compute the cross-correlation function, and it rec|»lres 
no attention from the user. In addition to CORR2, sub- 
routine SNRBPS requires GRAPH and DATAPT, which have al- 
ready been listed with the lowpass routine, and the FFT 
routine HARM. Actually, subroutines GRAPH and DATAPT are 
only used in either routine when a plot of the cross- 
correlation function is requested. Many FORTRAN compilers, 
however, require that the routines (or a dummy routine 
named GRAPH) be available in order to properly compile 
SNRBPS, which contains a calling statement for GRAPH. 
Similarly, CORR and C0RR2 are only executed in their res- 
pective routines if a delay estimate is to be found by the 

routines in forming the SNR estimate. 

The input parameter ISIGNR indicates whether a positive, 
negative, or an absolute value peak is to be found for the 
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cross-correlation function to select the j^eneral region 
where a search is performed to determine tjroup delay. The 
discussion given in Appendix A for SNRMSK explains the 
effect of ISIGNR on the search of the cross-correlation 
function. However, system gain is always considered to be 
positive in the bandpass case. The phase delay may be 
changed by 180° to provide sign changes. Therefore, in 
normal use, the user should always specify ISIGNR » 0. 
indicating the absolute value peak is to be found. 

iB some cases, signals may be analya^d where It Is 


known that theoretically equal magnitude positive and nega- 
tive cross-correlation peaks exist and that either peak 
would lead to a theoretically correct group delay estimate. 

A simple example of this is the sinusoidal modulation of a 
sinusoid used as a test signal for SNRBPS in the section on 
application to analog systems. It is preferable to use the 
smallest theoretically correct delay value available, since 
this procedure yields estimates based on the largest num- 
ber of samples. If the user knows a positive cross-correla- 
tion peak will occur for the smallest theoretically correct 
delay, it is helpful to specify ISIGNR - 1 to increase the 
likelihood of selecting this delay. Similarly, ISIGNR = 

-1 can be used if a negative peak is known to occur for the 
smallest theoretically correct delay. This procedure is 
occasionally useful, but it is entirely optional. The 
user may specify ISIGNR = 0 in all cases and obtian valid 
results. The option simply allows the user to obtain more 
reliable estimates in a few special cases. 
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APPENDIX D 

NOISE-EQUIVALENT BANDWIDTHS FOR THE DIGITAL FILTER 

For the purpose of computer simulation of conunun 1 - 
cation systems, a ^tt'nt'ral lowpass and bandpass bilimsir- 
Z dit;ital filter was desinn<*d with variabh* order ami 
critical frequencies. For varit>us filtt'r orders and 
a range of values of the filter critical frequencies, 
noise-equivalent bandwidths of the filter were measured 
using a tine dcxnain technique. A comparison was made 
between the noise-equivalent bandwidth of the lowpass 
digital filter and that of the corresponding lowpass 
analog filter. 

The filter was designed by first determining the 
transfer function of the Butterworth analog filter with 
the desired critical frequencies. Then standard bilinear— 
Z techniques were applied to obtain the correspond- 
ing digital filter. The filter was scaled so that 
H„.y. the maximum magnitude of the filter transfer func- 
tion, was equal to unity. 

The method used to measure the noise-equivalent 
bandwidth of the filter is based upon Parseval’s theorem, 
which states that [74] 


lH(f )raf - 


h^(t)dt 


where H(f) is the transfer function of the filter and 
h(t) is the impulse response. Since the filter is 
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causal. h(t) = 0 for t < 0. 


and 


I^OO ^ 

I |H(f)|=df . I 

^oo _ 

0 

By definition of B„, the 


(62) 


[75] 


noise-equivalent bandwidth. 


= HUOlV = [ |H(f,|2 

WAX 0 2II-_. J_« ' 


"MAX 0 


2«iAX 


df 


Since |»s been scaled to unity. 

8^*5 I |H(f)|2<,f. 


Substitution Of (62) into (63) yields 


(63) 




h^(t)dt 


for the noise-equivalent bandwidth. 

To measure the noise-equivalent bandwidth, a pulse 
Of unit area was applied to the digital filter. Rec- 
tangular integration was used by the computer to obtain 

B ~ 1 ? 2 

“ 2,1 h (nT)T 


ND 


n=0 


where T denotes the time between samples. The Integer 

K most be large enough to include all significant terms 

The value of K reguired depends upon how rapidly the 
summation converges. 



206 


Noise-equivalent bandwidths were computed for both 
lowpass and bandpass versions having' up to 20 poles. 

A wide ran^e of critical frequencies were eov.M-ed . and 
u rather extensive table was compiled for use with I h.- 
computer simulations of <-omnuin i cat ion systems. The table 
is n*.it rt'pri>diu‘ed lu'r*‘. .s i iuh' it is t>n 1 y aiiinali 
this particular software package applied to the lUM 070. 

In order to study the relationship between the noise- 
equivalent bandwidths for the analog and digital filters, 
an analysis was performed for several orders of the low- 
pass filter. The qualitative conclusions are easily 
extended to include the bandpass filter. 


The noise-equivalent bandwidth of the analog lowpass 
Butterworth filter can easily be shown to be [763 

^^3dB 

®NA 2Nsin(ir/2N) 

where f^^o denotes the 3dB frequency of the filter, and 
3ciB 

N is the filter order. 

Using the expressions for and ^ comparison 

was made of the noise-equivalent bandwidths of the low- 
pass analog filter and the corresponding lowpass digital 
filter. For filter orders 1 through 9, tests were per- 
formed over the range 0.0015 < R 1 0.25. where R is the 
3 dB frequency of the filter normalized by the sampling 
frequency. The sampling frequency, fg. was chosen to 
be unity for convenience. Thus, the normalized 3dB 
frequency. R. becomes equal to f ^3dB directly 
applicable to both digital and analog filters. 


I 


i 
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Results are given in Table VI for filter orders 
1, 2, 4, and 8. Curves are shown in Figure 34. It 
can be seen that the noise-equivalent bandwidths of 
the digital filter are less than thost? of the correspond- 
ing analog filters throughout this rang(?. now(>v<'r, th<‘ 
difference is small for small values ol R. The d i f f('rt*nc;<* 
becomes significant for large values of R and is clearly 
more pronounced for filters of low order . 

This behavior can be explained by remembering that 
the frequency response of the digital filter approaches 
zero as the frequency approaches fg/2. For frequencies 
which are very small with respect fo f^, the response 
curves of the analog and digital filters are nearly iden- 
tical. For frequencies near one half the sampling 
frequency, however, the digital filter frequency response 
begins to diverge downward from the analog response. If 
f„.„ is chosen so that R is small, the frequency responses 
of both filters are very low over the range where this 
divergence occurs, and the effect is small. Hence, 

Bjjjj is approximately equal to for very small values 

of R. If fojD is chosen so that R is large, the fre- 
quency response of the digital filter for frequencies 
above f^jo rolls off very rapidly. Thus, it immediately 

OClD 

diverges from the analog frequency response, which rolls 
off at the rate of 20N dB/ decade. It follows that ®ND 
and B„. are significantly different for large values of 

W A 

R. A filter of high order has a steeper frequency 
response curve than a low order filter. Therefore, a 
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table VI 
COMPARISON OF 

A,. DIGITAL NOISE-EQG.VALEHT BANn«IDTUS 


0.0015 
0.0050 
0.0150 
0.0500 
0.1500 
1 0.2500 



^^ND'^^^NA 


N=l 

N=2 

N=4 

0.995 

1.000 

1. 000 

0.985 

1.000 

1.000 

0.956 

0.999 

1.000 

0.870 

0.988 

1 0.998 

0.716 

0.939 

0.987 

0.637 

0.900 

0.975 


unstable 

unstable 

1.000 

1.000 

0.997 

0.994 


I , 

tn the nearest thousandth, 
All results are rounded to the 
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Figure 34* Comparison of Noise-Equivalent Bandwidths of Analog and Digital Lowpass 
Filters 
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greater value of R is required to show the difference 
between the values of and for a filter of high 

order. Exactly these effects can be seen in the plots. 

For values of R greater than about 0.25, accurate 
results for cannot be obtained using the techniques 

described because the integration unavoidably becomes 
coarse. Fixing the sampling frequency fixes T. which 
is the step size for the integration. For values of 
R between 0.25 and 0.5, the response, h(t), becomes 
large for t slightly greater than zero and decreases very 
rapidly to terms of negligible magnitude. Thus, the 
sum essentially converges after only a few terms, and 
the rectangular integration becomes inaccurate. This 
is not usually a problem, since the digital filter is 
most often designed to have critical frequencies less 

than one-fourth the sampling frequency. 

It is apparent that if an accurate value is re- 
quired for the noise-equivalent bandwidth of a digital 
filter, it is usually necessary to measure that band- 
width. Only for high order and high sampling rates do 
the well known analog noise-equivalent bandwidths yield 
accurate approximations tor the digital filter. Hence 
in the simulations, whenever a filter noise-equivalent 
bandwidth was required, the tabulated measurements were 
either utilized directly or used indirectly to Interpolate 
to the required values of the filter parameters. 
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APPENDIX E I 

DATA FOR COMMUNICATION SYSTEM SIMUEATIONS 

This appi'ndix contains tabU's of data cv' I I t'ot od < 

in performing simulations I'f both analog and digital i 

oommunioat ion systems. j 

Tables VII through XII contain data fi'^r the various j 

I 

cases of the FM communicat ion system siirulat ions . Each 
table is labeled with the modulation type, the value 
of the deviation ratio and the PLL natural frequency- 
used, as well as whether a hard-^limiter was utilized. 

Tables XIII through XVII contain the data which 
were obtained from the various digital ccmnunication 

system simulations. In each table, the particular system * 

which generated the data is given at the top. followed 

by a group of parameters for that simulation. Data | 

I 

concerning the SNR and the error probabi 1 i t ies appear • 

s 

last in each case. ■ 
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♦Flagged by SNRMSE for nonstationarity 

NOTE: All points except runs 10 and 14 were plotted, whether flagged or not. 
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TABLE XIII 

COHERENT ASK SIMULATION 


Carrier Frequency ^ 

Time-Bandwidth Product (4-Pole Filter) 

Sampling Frequency 

Number of Symbols Simulated 

Number of Samples Processed 

Number of Samples Used to Obtain SFR 


1 MHz 
2.244 
10 MHz 
1022 
30,700 
20,000 


Run 

No. 

SNR 

(dB) 

Pg(SNR) 

Pg( count) 

Number 

of 

Errors 

Phase Error 
Var . (Rad^) 


-4.21 

0.1781 

0.1703 

174 

0.0054 


-0.69 

0.0833 

0.0675 

69 

0.0082 


1.37 

0.0397 

0.0352 

36 

0.0122 


4.17 

0.00774 

0.00783 

8 

0.0171 


4.97 

0.00398 

0.00196 

2 

0.0172 

6 

-3.15 

0.1485 

0.1184 

121 

0 . 00910 

7 

0.46 

0.0572 

0.0470 

48 

0.0196 

A 

2.61 

0.0216 

0.0166 

17 

0 . 0249 

w 

5.60 

0.00217 

— 

0 

0,0219 

m 

6.48 

0.000797 


0 

0 . 0205 



I 



Runs 1-5 utilize filtered signal-plus-noise 
Runs 6-10 add filtered noise to the signal 
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TABLE XIV 


CUHERENT FSK SIMULATION 

Frequency for Symbol 1 
Frequency for Symbol 2 

Product (4-Pole Filter) 
Sampling Frequency 

Number of Symbols Simulated 
Number of Samples Processed 
Number of Samples Used to Obtain SNR 


1.00 MHz 
1.67 MHz 
4.458 

15 MHz 

1023 

46,100 

20.000 


Run 

No. 

SNR 

(dB) 

yv y\ 

Pj(SNR) 

Pg( count) 

Number 

of 

Errors 

Phase Error 
Var. (Rad^) 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

-7.20 

-3.60 

-1.47 

1.50 
2.37 

1 -7.41 

i -3.71 
-1 . 52 

1.51 
2.40 

0.1782 

0.0816 

0.0373 

0.00605 

0.00277 

0.1842 

0.0842 

0.0382 

0.00600 

0.00269 

0.1750 
0.0763 
0.0381 
0.00196 
0.000978 
j 0.1779 
0.0841 
0.0313 
0.00599 
0.00269 

179 

78 

39 

2 

1 

182 

86 

32 

6 

1 

0.033 

0.015 

0.010 

0.006 

0.005 

0.036 

0.015 

0.009 

0.004 

0.004 

0.030 
0.015 
0.011 
0.008 
0.008 
0.036 
0.015 
! 0.009 
0.004 
0.004 


NOTE: 


Runs 1-5 utilize filtered signal-plus-noise 
Runs 6-10 add filtered noise to the signal 
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TABLE XV 


NONCOHERENT FSK SIMULATION 
WITHOUT RAYLEIGH FADING 


Frequency for Symbol 1 
Frequency for Symbol 2 
Time-Bandwidth Product (Predetection 
4-Pole Filter) 

Saiq>ling Frequency 
NuBber of Symbols Simulated 
Number of Siuq>les Processed 
Number of SaBq>les Used to Obtain SNR 


1.00 MHz 
1.67 MHz 

6.390 
15 MHz 
1023 
46.100 

20.000 


Run 

No. 

■■ 


Pg( count) 

Number 

of 

Errors 

1 

-5.76 

0.2141 

0.2180 

223 

2 

-2.16 

0.0717 

0.0850 

87 

3 

-0.02 

0.0208 

0.0303 

31 

4 

0.94 

0.0095 

0.0117 

12 
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TABLE XVI 

NONCOHERENT FSK SIMULATION IN 
SINGLE RAYLEIGH FADING CHANNEL 

All system and simulation parameters are identical to 
those given in Table XV 


Run 

No. 



Pg( count ) 

Number 

Er?irs 

1 

-8.00 

0.3319 

0.3509 

359 

2 

-2.98 

0.1917 

0.2190 

224 

3 

2.02 

0.0822 

0.0958 

98 

4 

4.06 

0.0547 

0.0499 

51 

5 

6.97 

0.0296 

0.0352 

36 

6 

11.79 

0.0102 

0.0978 

10 


TABLE XVII 

DIVERSITY TRANSMISSION IN 
THE RAYLEIGH FADING CHANNEL 


All system and simulation parameters are identical to 

those given in Table XV. Channel noise was adjusted to 

maintain (E /N ) = 15 dB 
s' o 


Run 

No. 


Order of 
Diversity 

SNR (dB) 
Per 

Channel 

Theo.P„ 
.. E 

at 

15 dB 

P„(count ) 

Ei 

Number 

of 

Errors 

1 

6.97 

0.0297 

0.0323 

33 

2 

4.06 

0.00910 

0.00978 

10 

3 

2.89 

0.00448 

0.00489 

5 

4 

1.83 

0.00283 

0 . 00489 

5 
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APPENDIX F 

THE EFFECT OF SNR ESTIMATION ERRORS ON P^ ESTIMATES 

E 

When the SNR (signal-to-noise ratio) is estimated at 
a point in a particular simulation, errors in the estimate 
generally occur. The error, e, in the SNR estimate, 

SNR, can be defined by 

/V 

SNR = (l+e)(SNR) (64) 

so that e=0 corresponds to an exactly correct estimate. 

Many times in applications to digital systems, the prob- 

ability of symbol error, P^, is estimated using SNR. Thus, 
the error in SNR produces an error in the symbol error 

/N 

probability estimate, P^, and it is useful to evaluate the 

^ /s, 

error in Pg resulting from the error in SNR. The sensitiv- 
ity of Pg to the measurement error, c, is dependent upon 

/V ys 

the mapping from SNR to Pj, for the particular system in- 
vo 1 ved . 

A number of coherent systems operating in additive 
Gaussian noise environments can be studied using a general 
mapping of the form 

Pr = f erfc(/z) (65) 

where 

“ 2 

erfc (X) = ^ [ e'“ da. 

/if ■' X 


The quantity z represents the product of a system depen- 
dent constant, 6, with the ratio of the average energy per 
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transmitted symbol, E,. to the noise power spectral density 

N . Therefore, (65) may be written 
o 


It was shown In the section on application to dipital 
systems that 


^ W (SNRKV^) 


where (T^B„) is the time-bandwidth product for the system. 
Therefore the general mapping may be written 

p„ e i erfc /6<SHR)(T • 

E ^ 

This relation may be applied to many communication systen.. 
If B = ^ la baed, the mapping is valid for a coherent bi- 
nary ASK or FSK system. With 6 - 1. the mapping applies 
to a phase- reversal keying system [77]. The mapping can 
also be applied to a quadriphase phase-shift keying sys- 


tem [78]- 

Since many systems may be represented by this general 
mapping, it Is helpful to Investigate the effect of errors 
in s5b on the resulting error probability estimate obtained 


using 


p = i erfc /e(SNR)(TgBj^). 
£ Z 


Substituting from (64) yields 
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P., = 


_ 1 


2 erfc ./6(SNR)(T B^)(l+e). 

o n 


Noting 


E 

2 = B jp = S(SNR)(T B^) 
o ® ” 


yields 

p£ “ 2 ®rfc/z( 1+e ) (66) 

which gives a general form for analysis. Plots of (66) 
are given in Figure 35 for values of e from -0.5 to 0.5. 
The curve with e=0 represents the mapping when no error 
is present. The distance along the ordinate between this 
curve and any other curve is the error in P^. which corres- 
ponds to the value of e specified. It can be seen that 
for large vlaues of z, the curves diverge, and the slopes 
of the curves increase. This leads to larger errors in 

/V 

^ given value of e as z is increased. Figure 35 
IS not particularly convenient for studying the error in 

i(N 

Pg, so another representation is helpful. 

In many applications P^. estimates are typically 
specified by giving the order of magnitude. For example 
Pg = 0.871x10 would usually be simply referred to as 
Pg=10 . Often the order of magnitude is an adequate 

specification for P^. and for specifying error in P_. Per- 
centage error measurements are not convenient to use for 
specifying the error in P^ in most cases. An appropriate 
way to specify the error in Pj, is provided by defining 


Estimated Probability of Symbol Error 







Using this expression. >-0 corresponds to 
Indicating zero error. Additionally. 
yields y = l.(Pg/Pj,) = 1/10 yields > = -1. and 
(P /P^) = l/lOO yields i = -2 . and so on. Thus. 

ill tl* 

Y is the order of magnitude error in the estimate, as 
desired. 

A plot of Y versus 10 ‘LOG^q(z) is giver in Figure 
36. This plot illustrates the sensitivity of the error 
in P-e to the error in the estimate for z. The divergence 
of the curves for increasing values of z indicates how a 
given value of e produces larger errors in P^. for larger 
values of SNR. The box shown in the figure encompasses 
the region where the coherent ASK and FSK simulations were 
performed. This region is shown in Figure 37. On this 
plot the results obtained in the coherent ASK and FSK 
simulations are shown. The direct error counts are used 
to define the probability of symbol error, Pj,. and the 
estimate for the symbol error probability computed using 
SNR estimation defines Pg. For each point where less 
than five errors were counted, the number of errors is 
shown next to the point. The data for the simulations are 
tabulated in Appendix E. 

It can be seen that excellent results were obtained. 
In all cases where more than five error.s were counted, so 
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